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.

Solución análítica de ecuaciones diferenciales ordinarias

En esta entrada usaremos MatLab para encontrar la solución analítica de una ecuación diferencial ordinaria. Para ello utilizaremos la función dsolve.
--------------------------------------------------------------------------------------------------------------------------------
dsolve
Soluciona ecuaciones diferenciales ordinarias y sistemas de ecuaciones diferenciales ordinarias
Sintaxis
S = dsolve(eqn)
S = dsolve(eqn,cond)
S = dsolve(eqn,cond,Name,Value)
Y = dsolve(eqns)
Y = dsolve(eqns,conds)
Y = dsolve(eqns,conds,Name,Value)
[y1,...,yN] = dsolve(eqns)
[y1,...,yN] = dsolve(eqns,conds)
[y1,...,yN] = dsolve(eqns,conds,Name,Value)
Descripción 
En esta entrada solo describiremos las dos primeras.
1). S = dsolve(eqn): Soluciona la ecuación diferencial ordiaria eqn. Aquí eqn es una ecuación simbólica que contiene diff para indicar derivadas.
 2). S = dsolve(eqn,cond): Soluciona la ecuación diferencial ordinaria eqn con la condición inicial o condición de contorno cond.

Las variables simbólicas se pueden definir con la función syms. Alternativamente, se puede usar un string con la letra D para indicar derivadas.
----------------------------------------------------------------------------------------------------------------------------- 
Ejemplo 1: Obtengamos la solución particular de la ecuación diferencial ordinaria de primer orden:

                            y'= sin(2*x)-y*tan(x);   y' = dy/dx,
con la condición inicial
                            y(0)=1.
----------------INICIO------------------------------------------------------------------------------------------------------

%---creamos las variables simbólicas con syms:
syms y(x)
%----Construimos la ecuación diferencial y la colocamos en la función dsolve:
y(x) = dsolve(diff(y,x) == sin(2*x)-y*tan(x),y(0)==1);
%----Evaluamos y para valores de x entre 0 y 4 a pasos de 0.1 y los valores lo asignamos a la variable y1: 
x = 0:0.1:4;
y1=eval(y);
%----graficamos y1 usando la función plot:
figure
plot(x,y1)
grid on
%---Evaluamos y para el valor x = 1 usando eval: 
 x=1;
yo=eval(y);
%---graficamos el punto  (x,yo) sobre la gráfica de y(x), de color rojo 'r' y marca 'o' :
hold on
plot(x,yo,'ro')
xlabel('x')
ylabel('y')
grid on

-------------------FIN-----------------------------------------------------------------------------------------------------------

Los resultados se muestran a continuación:
La solución es:
y(x) = 3*cos(x)-cos(2x)-1
Evaluada para x=0 resulta:
yo = 1.0371

La gráfica resultante es:


------------------------------------------------------------------------------------------------------------------------------
Ejemplo 2: Obtengamos la solución particular de la ecuación diferencial ordinaria de segundo orden:
                          h'' + 2h'+2h = 0;  h' = dh/dt,
con la condición inicial
                           h(0)=0  ;  h''(0)=15 .
----------INICIO-------------------------------------------------------------------------------------------------------
%---creamos las variables simbólicas con syms:
syms h(t)
%----Construimos la ecuación diferencial y la colocamos en la función dsolve. Usamos D en lugar de diff para definir la ODE colocándola como string (entre comillas simples):
h(t) = dsolve('D2h == -2*Dh-2*h', h(0)==0, 'Dh(0)==15')
%----Evaluamos h para 100 valores de t entre 0 y 10 y lo asignamos a la variable h1:
N=100;
t = linspace(0,10,N);
h1=eval(h);
%----graficamos h1 usando la función plot:
figure
plot(t,h1)
grid on
%---Evaluamos h para el valor t = 2:
t=2;
ho=eval(h);
%---graficamos el punto (0,ho) sobre la gráfica de h(t), de color rojo 'r' y marca '*' :
hold on
plot(t,ho,'r*')
xlabel('t')
ylabel('h')
grid on
-------------FIN-------------------------------------------------------------------------------------------------------------

Los resultados se muestran a continuación:
La solución es:
h(t) = 15*exp(-t)*sin(t)
Evaluada para t=2 resulta:
ho = 1.8459
La grafica resultante es:

Gráficas de campos direccionales con MatLab

Esta rutina en Matlab realiza la gráfica de un campo direccional a partir de la ecuación diferencial :

y'=f(x,y)=-5*x.^4.*y.^2.

Sobre el campo de direcciones se realiza la gráfica de la solución particular:

y = 1./(x.^5+1).

%=====INICIO========
%-1)---Región--------
x1=0;
x2=2;
y1=0;
y2=2;
N = 20;
x = linspace(x1,x2,N);
y = linspace(y1,y2,N);
[X,Y]= meshgrid(x,y);

%-2)---calculo del angulo de inclinación de cada vector
m = -5*X.^4.*Y.^2;
T = atan(m);

%-3)--componentes de cada vector
U = cos(T);
V = sin(T);

%-4)--Dibujo de los vectores tangentes
h = quiver(X,Y,U,V);

%-5)--superposición en el campo de direcciones de la solución particular
hold on
yp = 1./(x.^5+1);
plot(x,yp)
grid on
%===FIN========

El resultado es el siguiente:

Fig.1. Campo de direcciones (azul). Solución particular (rojo)