Ayuda con resolucion de ecuacion diferencial
Publicado por
Mario Garcia (1 intervención) el 04/10/2011 05:25:21
Gracias por tu ayuda, la verdad que la ecuación, cuando la resolvía en la computadora de la universidad me salia otra respuesta a la que tu me estas dando, pero cuando he intentado resolver esta ecuación en casa y evaluarla, me ha salido el mismo resultado y no he tenido ningún problema.
Justo esa ecuación diferencial la he aproximado por los métodos numéricos de Euler, Euler Modificado, Runge-Kutta-2 y runge-Kutta-4. Si no es mucha molestia ¿ crees que le podrías dar un vistazo a mi código y opinar?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Problema 1:
%1a)
%La ecuacion a resolver por los Metodos de Euler y Euler modificado
%es: y''-1/4*y=1/4*t*e^(t/2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Metodo de Euler:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%y'' -1/4*y=1/4te^t/2 se convierte a un sitema de primer orden de dos
%variables. para esto definimos y'=x ^ x'=y'', de esta manera podemos armar
%el sistema para x' e y' :
ymetodo=@(y,x,t) x;
xmetodo=@(y,x,t) (1/4).*y+(1/4).*t.*exp(t/2);
y1(1)=1;
x1(1)=0;
h=0.1;
tend=10;
t=0:h:10;
nstep=tend/h;
%se resuelve el sistema a través de una resolución iterativa:
for i=1:nstep;
y1(i+1)=y1(i)+h*ymetodo(y1(i),x1(i),t(i));
x1(i+1)=x1(i)+h*xmetodo(y1(i),x1(i),t(i));
end;
%Con esto se ha terminado de calcular el la trayectoria de Y aproximada con
%el método de Euler.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Euler Modificado:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ym(1)=y1(1);
xm(1)=x1(1);
%Proceso de Iteración
for a=1:nstep;
t1=ymetodo(xm(a),ym(a),t(a));
t2=xmetodo(xm(a),ym(a),t(a));
ym(a+1)=ym(a)+h*ymetodo(ym(a)+(h/2)*t1,xm(a)+(h/2)*t2,t(a)+h/2);
xm(a+1)=xm(a)+h*xmetodo(ym(a)+(h/2)*t1,xm(a)+(h/2)*t2,t(a)+h/2);
end;
figure('name','Gráfico: Aproximacion Euler y Euler Modificado');
subplot(1,1,1);plot(t,y1,t,ym);legend('Euler','Euler-Modificado');axis tight;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%1b)
%resolver y graficar el mismo sistema por los metodos de RK-2 y RK-4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%RK-2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yr2(1)=y1(1);
xr2(1)=x1(1);
for b=1:nstep
yr2(b+1)=yr2(b)+(h/2)*(ymetodo(yr2(b),xr2(b),t(b))+ymetodo(y1(b+1),x1(b+1),t(b+1)));
xr2(b+1)=xr2(b)+(h/2)*(xmetodo(yr2(b),xr2(b),t(b))+xmetodo(y1(b+1),x1(b+1),t(b+1)));
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%RK-4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yr4(1)=y1(1);
xr4(1)=x1(1);
for c=1:nstep
temp1y=ymetodo(yr4(c),xr4(c),t(c));
temp1x=xmetodo(yr4(c),xr4(c),t(c));
temp2y=ymetodo(yr4(c)+(h/2)*temp1y,xr4(c)+(h/2)*temp1x,t(c)+h/2);
temp2x=xmetodo(yr4(c)+(h/2)*temp1y,xr4(c)+(h/2)*temp1x,t(c)+h/2);
temp3y=ymetodo(yr4(c)+(h/2)*temp2y,xr4(c)+(h/2)*temp2x,t(c)+h/2);
temp3x=xmetodo(yr4(c)+(h/2)*temp2y,xr4(c)+(h/2)*temp2x,t(c)+h/2);
temp4y=ymetodo(yr4(c)+h*temp3y,xr4(c)+h*temp3x,t(c)+h);
temp4x=xmetodo(yr4(c)+h*temp3y,xr4(c)+h*temp3x,t(c)+h);
yr4(c+1)=yr4(c)+(h/6)*(temp1y+2*temp2y+2*temp3y+temp4y);
xr4(c+1)=xr4(c)+(h/6)*(temp1x+2*temp2x+2*temp3x+temp4x);
end;
figure('name','Gráfico: Aproximacion RK-2 y RK-4');
subplot(1,1,1);plot(t,yr2,t,yr4);legend('RK-2','RK-4');axis tight
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%1c)
%Obtener la trayectoria verdadera y los errores de los cuatro metodos
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%caluclo de la solucion exacta
yv=dsolve('D2y-(1/4)*y=(1/4)*t*exp(t/2)','y(0)=1','Dy(0)=0');
yv=eval(yv);
%calculo del Error Euler
y1e=100.*abs(yv-y1)./yv;
%calculo del Error Euler modificado
yme=100.*abs(yv-ym)./yv;
%calculo del Error RK-2
yr2e=100.*abs(yv-yr2)./yv;
%calculo del Error RK-4
yr4e=100.*abs(yv-yr2)./yv;
%Graficamos la trayectoria Verdadera y las 4to Aproximaciones asi como los
%errores
figure('name','Gráficos:Trayectorias y Errores');
subplot(2,1,1);plot(t,yv,t,y1,t,ym,t,yr2,t,yr4);legend('Y-verdadero','Y-Euler','Y-Euler Modificado','Y-RK2','Y-RK4');axis tight;
subplot(2,1,2);plot(t,y1e,t,yme,t,yr2e,t,yr4e);legend('Euler-E','Euler Modificado-E','RK2-E','RK4-E');axis tight
% FIn de pregunta 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
