Diferencia entre revisiones de «Reacciones con autocatálisis 4-C»
De MateWiki
(→Resolución del P.V.I mediante los métodos numéricos: Euler, trapecio y Rounge-Kutta) |
(→Resolución del P.V.I mediante los métodos numéricos: Euler, trapecio y Rounge-Kutta) |
||
| Línea 77: | Línea 77: | ||
hold off | hold off | ||
}} | }} | ||
| + | <br /> | ||
| + | [[Image:JJ2.png|500px|thumb|left|figure 1]] [[Image:JJ2a.png|500px|thumb|center|Zoom figure 2]] | ||
| + | |||
| + | <br/> <br/> | ||
| + | |||
[[Archivo:JJ2.png|300px|miniaturadeimagen|centro|izquierda|figure2]] | [[Archivo:JJ2.png|300px|miniaturadeimagen|centro|izquierda|figure2]] | ||
Revisión del 14:16 4 mar 2015
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Visualizaci´on de campos escalares y vectoriales en elasticidad. Grupo 4-C |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2014-15 |
| Autores | Ignacio Mollá Carcaño, Jorge Javier Rodríguez Anzules, Claudia Jalón Manzano, Pablo Revuelta Aragón, Pedro Torrecilla Sánchez, Alejandro Martínez Gamonal |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
1 Introdución:Reacción bimolecular: A+B[math]\longrightarrow[/math] C
Estamos ante una reacción química irreversible autocatalítica, [math]A+B \longrightarrow C [/math]. Esta reacción obviamente cumplirá el principio de conservación de la masa, es decir, la suma de las masas de los reactivos deberá ser igual a la de los productos, o dicho de otra forma, la velocidad de destrucción es igual a la de creación. Además sabemos que la velocidad de reacción es igual al producto de los reactivos multiplicado por la constante cinética.
[math] X_{0}=1 mol/l[/math] e [math]Y_{0} =0.001 mol/l[/math] ,siendo [math]Y=B[/math] , [math]X=A [/math] y [math] K=1 [/math]
[math] y'(t)=kx(t)y(t)[/math] [math] y'(t)=y(t)(1.01-y(t)) [/math]
2 Resolución del P.V.I mediante los métodos numéricos: Euler, trapecio y Rounge-Kutta
clear all
clf
clc
%Trabajo 2: Reacciones con autocatálisis
%Datos
h=0.1;
t0=0;
y0=0.01;
x0=1;
tN=10;
%Calculamos los vectores(Apartados 2 y 3)
t=t0:h:tN; % vector de tiempos
N=(tN-t0)/h;
y=zeros(1,length(t)); %vector de concentración de B según euler
x=zeros(1,length(t)); %vector de concentración de A según euler
z=zeros(1,length(t)); %vector de concentración de B según trapecio
r=zeros(1,length(t)); %vector de concentración de B según rounge-kutta
e=zeros(1,length(t)); %vector de concentración de B exacto
y(1)=y0;
x(1)=1;
z(1)=y0;
r(1)=y0;
for i= 1:N;
%Euler
y(i+1)=y(i)+h*(y(i)*(1.01-y(i)));
x(i)=1.01-y(i);
%Trapecio
z(i+1)=((1.01*(h/2)-1)+sqrt((1-1.01*(h/2))^2-2*h*((h/2)*(z(i)^2)-(1+1.01*(h/2))*z(i))))/h;
%Rounge-kutta
K1=r(i)*(1.01-r(i));
K2=(r(i)+1/2*K1*h)*(1.01-(r(i)+1/2*K1*h));
K3=(r(i)+1/2*K2*h)*(1.01-(r(i)+1/2*K2*h));
K4=(r(i)+(K3*h))*(1.01-r(i)+(K3*h));
r(i+1)=r(i)+(h/6)*(K1+2*K2+2*K3+K4);
%Exacta
e(i+1)=(1.01/(100*exp(-1.01*t(i+1))+1));
end
figure(1)
hold on
plot(t,y);
plot(t,x,'r');
xlabel('tiempo(s)'); ylabel('concentración(mol/l)');
legend('Elemento B','Elemento A','Location','best');
title('Gráfico de las concentraciones de A y B');
hold off
figure(2)
hold on
plot(t,y,'b') %euler
plot(t,z,'g') %trapecio
plot(t,r,'k') %rounge-kutta
plot(t,e,'r') %exacta
xlabel('tiempo(s)'); ylabel('concentración de B(mol/l)');
legend('Euler','Trapecio','Rounge-kutta','Exacta','Location','best');
title('Gráfico de la concentración de B según los diferentes métodos')
hold off
figure(3)
%gráfico de errores
ey=zeros(1,length(t));
et=zeros(1,length(t));
er=zeros(1,length(t));
ey(1)=0; et(1)=0; er(1)=0;
for i=1:N
ey(i+1)=abs(y(i+1)-e(i+1)); %errores euler
et(i+1)=abs(z(i+1)-e(i+1)); %errores trapecio
er(i+1)=abs(r(i+1)-e(i+1)); %errores rounge-kutta
end
hold on
plot(t,ey,'b')
plot(t,et,'g')
plot(t,er,'k')
legend('Errores euler','Errores trapecio','Errores rounge-kutta','Location','best');
xlabel('tiempo(s)'); ylabel('error cometido(mol/s)');
title('Grafico de errores')
hold off%Por sistemas de ecuaciones (Apartado 4)
xx=zeros(1,length(t)); %vector de concentración de A según euler
yy=zeros(1,length(t)); %vector de concentración de B según euler
ra=zeros(1,length(t)); %vector de concentración de A según rounge-kutta
rb=zeros(1,length(t)); %vector de concentración de B según rounge-kutta
xx(1)=x0;
yy(1)=y0;
ra(1)=x0;
rb(1)=y0;
for i=1:N
%Sistema con Euler
xx(i+1)=xx(i)-h.*(xx(i).*yy(i));
yy(i+1)=yy(i)+h.*(xx(i).*yy(i));
%Sistema con Rounge-kutta de orden 4
K1a=-1*ra(i)*rb(i);
K1b=1*ra(i)*rb(i);
K2a=-1*(ra(i)+1/2*K1a*h)*(rb(i));
K2b=1*(rb(i)+1/2*K1b*h)*(ra(i));
K3a=-1*(ra(i)+1/2*K2a*h)*(rb(i));
K3b=1*(rb(i)+1/2*K2b*h)*(ra(i));
K4a=-1*(ra(i)+K3a*h)*(rb(i));
K4b=1*(rb(i)+K3b*h)*(ra(i));
ra(i+1)=ra(i)+(h/6)*(K1a+2*K2a+2*K3a+K4a);
rb(i+1)=rb(i)+(h/6)*(K1b+2*K2b+2*K3b+K4b);
end
figure (4)
hold on
plot(t,xx,'b')
plot(t,yy,'r')
plot(t,ra,'g')
plot(t,r,'k')
xlabel('tiempo(s)'); ylabel('concentración(mol/l)');
legend('Concentracion de A con Euler','Concentracion de B con Euler','Concentración de A con Rounge-kutta','Concentración de B con Rounge-kutta','Location','best');
hold off
