Métodos para ecuaciones diferenciales de primer orden: Metodo de Euler, Euler mejorado Y Runge-Kutta 4° orden

Con esta entrada comenzaremos con el uso de métodos numéricos para resolver ecuaciones diferenciales. Comenzaremos con los sencillos métodos de Euler y Euler mejorado para la solución de ecuaciones diferenciales de primer orden con valor inicial y finalizaremos con el método de Runge-Kutta de cuarto orden.

 Ecuación diferencial de primer orden con valor inicial:

 y' = f(x,y),         y(xo) = yo.

Utilizaremos los métodos de Euler para resolver la ecuación diferencial
   
y' = y + x,      y(0) = 0.

Esta ecuación tiene solución analítica    y  =  -1 + exp(x) - x.

A modo de comparación graficaremos la solución encontrada con el método numérico y la solución analítica.
 --------------------------------------
1). METODO DE EULER
 ---------------------------------------

yn+1 = yn + hf(xn, yn);              n = 0,1,2,3...,N 

El programa realizado en MatLab es el siguiente:

%-----INICIO-----
clc
close all
clear all
h = 0.2; %---paso
N=10;    %---numero de valores a calcular
%----valores iniciales
xo = 0;
yo = 0;
%---Definición del vector "y" de longitud N e inicialmente con ceros.
y = zeros(1,N);
y(1)=yo; %--valor inicial asignado al primer elemento de y.
yn = yo;
for n=1:N-1;
    xn = xo+(n-1)*h;
    f = (yn+xn); %--funcion f que depende de la ODE a resolver
    y(n+1) = yn + h*f;
    yn = y(n+1);
end
%---graficas
x = 0:h:h*(N-1);
ya = -1 + exp(x)-x; %--solución analítica
plot(x,ya,'k',x,y,'r')
hold on
plot(x,ya,'k*',x,y,'ro')
xlabel('x')
ylabel('y')
legend('analytical solution','Euler Method')
grid on
%---FIN

Resultado

 
 ----------------------------------------------------------
2). METODO DE EULER MEJORADO
 ----------------------------------------------------------
 y*n+1 = yn + hf(xn + yn)

yn+1 = yn + (1/2)h[f(xn, yn) + f(xn+1 , y*n+1)];              n = 0,1,2,3...,N 

El programa realizado en MatLab es el siguiente:

%-----INICIO-----
clc
close all
clear all
h = 0.2;
N=10;  
xo = 0;
yo = 0;
y = zeros(1,N);
y(1)=yo;
y1 = yo;
for m=1:N-1;
    x1 = xo+(m-1)*h;
    x2 =x1 +h;
    f1 = y1+x1;
    y2 = y1+h*f1;
    f2 = y2+x2;
    y(m+1) = y1 + 0.5*h*(f1+f2);
    y1 = y(m+1);
end

 x = 0:h:h*(N-1);
ya = -1 + exp(x)-x;

figure
plot(x,ya,'k',x,y,'r')
hold on
plot(x,ya,x,y,'*r')
legend('Analytical solution', 'Improved Euler Method')
xlabel('x')
ylabel('y')
grid on
%---END

Resultado



 ---------------------------------------------------
2). MÉTODO DE RUNGE-KUTTA
 ---------------------------------------------------

 k1 = hf(xn, yn)
 k2 = hf(xn + 0.5h, yn + 0.5k1)
 k3 = hf(xn + 0.5h, yn + 0.5k2)
 k4 = hf(xn + h, yn + k3)
 xn+1 = xn+h
 yn+1 = yn + (1/6)(k1 + 2k2 + 2k3 + k4)
 n = 0,1,2,3...,N 

El programa realizado en MatLab es el siguiente:

yr = zeros(1,N);
yr(1)=yo;
y1 = yo;
for r=1:N-1;
    x1 = xo+(r-1)*h;
    k1 = h*(x1+y1);
    k2 = h*(x1+0.5*h+y1+0.5*k1);
    k3 = h*(x1+0.5*h+y1+0.5*k2);
    k4 = h*(x1+h+y1+k3);
    yr(r+1) = y1 + (1/6)*(k1+2*k2+2*k3+k4);
    y1 = yr(r+1);
end 
x = 0:h:h*(N-1);
ya = -1 + exp(x)-x;
figure
plot(x,ya,'k',x,yr,'r')
hold on
plot(x,ya,'k.',x,yr,'ro')
legend('Analytical solution', 'Runge-Kutta Method')
grid on
xlabel('x')
ylabel('y')
Resultado


Comparación gráfica:  Solución analítica, método de Euler, Método de Euler mejorado Y método de Runge-Kutta.