Solución numérica de ecuaciones diferenciales (I)

Método de Euler

Vamos aresolver la ecuación diferencial de primer orden

d x d t = f ( t , x )

con con la condición inicial de que en el instante t0 la posición es x0

La primera derivada nos permite conocer la posición xi+1 en el instante ti+1, a partir de la posición xi en el instante ti de acuerdo a la fórmula siguiente. La línea de color rojo es la tangente a la curva en el instante ti

xi+1=xi+f(ti,xi)h

El procedimiento de Euler produce un error que se acumula a cada paso h de integración, que es el segmento en color azul que une los dos puntos en la figura.

Escribimos una función denominada euler, a la que le pasaremos:

y nos devolverá un vector t y su correspondiente vector x.

function [t,x] =euler(f,t0,tf,x0,n)
      h=(tf-t0)/n;
      t=t0:h:tf;
      x=zeros(n+1,1); %reserva memoria para n+1 elementos del vector x
      x(1)=x0;
      for i=1:n
          x(i+1)=x(i)+f(t(i),x(i))*h;
      end
end

Supongamos que queremos integrar la ecuación diferencial

dx dt =cost

con las condición inicial t=0, x=0.

x0= 0 t cost·dt x=sint

Tomamos un intervalo h=π/6, y construimos la siguiente tabla

t

dx dt =cost

x(Euler)

x=sin t

0 1 0 0
π/6 0.866 0.523 0.5
π/3 0.5 0.977 0.866
π/2 0 1.239 1
2π/3 -0.5 1.239 0.866
5π/6 -0.866 0.977 0.5
π 0.523 0

Para aplicar el método de Euler precisamos de un paso h pequeño, incluso así los errores se van acumulando y al cabo de cierto tiempo la diferencia entre el valor exacto y el calculado es grande.

Escribimos un script en el que definiremos la función f(t,x), las condiciones iniciales y llamaremos a la función euler. Finalmente, representaremos gráficamente la solución exacta y la obtenida aplicando el método de Euler

tf=input('tiempo final, tf: ');
n=input('número de pasos, n: ');
f=@(t,x) cos(t); 
%condiciones iniciales
t0=0; 
x0=0;
[t,x]=euler(f,t0,tf,x0,n); 

hold on
plot(t,x,'b')
y=sin(t);
plot(t,y,'r')
xlabel('t')
ylabel('x')
grid on
legend('aproximada','exacta')
title('dx/dt=cost')
hold off

En la ventana de comandos corremos el script

tiempo final, tf: pi
número de pasos, n: 40

Apreciamos diferencias entre la solución exacta y la obtenida mediante integración numérica por el método de Euler

Método de Runge-Kutta

En esta sección vamos a estudiar la aplicación del método de Runge-Kutta a:

Ecuación diferencial de primer orden

Sea una ecuación diferencial de primer orden, con la condición inicial de que en el instante t0 el valor inicial de x es x0

Se elige una anchura de paso h y se calculan cuatro números k1, k2, k3, k4 de acuerdo con el procedimiento esquematizado en la tabla adjunta. Según el procedimiento ordinario de Runge-Kutta, a partir del valor de x en el instante t se determina el valor de x en el instante t+h mediante la fórmula que figura en la última fila de dicha tabla.

dx dt =f(t,x)

k 1 =h·f(t,x) k 2 =h·f( t+ 1 2 h,x+ 1 2 k 1 ) k 3 =h·f( t+ 1 2 h,x+ 1 2 k 2 ) k 4 =h·f( t+h,x+ k 3 )

x(t+h)=x(t)+ 1 6 ( k 1 +2 k 2 +2 k 3 + k 4 )

Definimos la función rk_1 que resuelve la ecuación diferencial de primer orden, cuando le pasamos:

y nos devolverá un vector t y su correspondiente vector x.

function [t,x] =rk_1(f,t0,tf,x0,n)
    h=(tf-t0)/n;
    t=t0:h:tf;
    x=zeros(n+1,1); %reserva memoria para n elementos del vector x
    x(1)=x0;
    for i=1:n
        k1=h*f(t(i),x(i));
        k2=h*f(t(i)+h/2,x(i)+k1/2);
        k3=h*f(t(i)+h/2,x(i)+k2/2);
        k4=h*f(t(i)+h,x(i)+k3);
        x(i+1)=x(i)+(k1+2*k2+2*k3+k4)/6;
     end
end

Considérese el circuito en serie de la figura. Inicialmente el condensador está descargado. Si se cierra el interruptor I la carga empieza a fluir produciendo corriente en el circuito, el condensador se empieza a cargar. Una vez que el condensador adquiere la carga máxima, la corriente cesa en el circuito.

R dq dt = V ε q C 0 q dq C V ε q = 1 RC 0 t dt q=C V ε ( 1exp( t RC ) )

Escribimos un script para que realice las siguientes tareas:

  1. Establezca, mediante comandos input:
  2. Fije las condiciones iniciales, en el instante inicial t=0, el condensador está descargado x=0.
  3. Defina la función f(t,x),
  4. Llame al procedimiento numérico rk_1
  5. Mediante el comando plot realice una representación gráfica de la solución numérica
  6. Realice una representación gráfica de la solución exacta

Ejemplo: R=2.0, C=0.8, y tf=10.

V0=10;
R=input('Resistencia R: ');
C=input('Capacidad C: ');
tf=input('tiempo final, tf: ');
n=input('número de pasos, n: ');

f=@(t,x) V0/R-x/(R*C); 
%condiciones iniciales
t0=0; x0=0; 
[t,x]=rk_1(f,t0,tf,x0,n); 

hold on
plot(t,x,'b')
y=C*V0*(1-exp(-t/(R*C)));
plot(t,y,'r')
xlabel('t')
ylabel('q')
grid on
legend('aproximada','exacta','Location','Southeast')
title('Carga del condensador')
hold off

En la ventana de comandos corremos el script

Resistencia R: 2
Capacidad C: 0.8
tiempo final, tf: 10
número de pasos, n: 50

rk_1

No se aprecia diferencia entre la solución exacta y la numérica, aplicando el procedimiento de Runge_Kutta.

Sistema de dos ecuaciones diferenciales de primer orden

El procedimiento de Runge-Kutta es igualmente efectivo en la resolución de un sistema de dos ecuaciones diferenciales de primer orden.

dx dt =f( t, x, y) dy dt =g( t, x, y)

El procedimiento de aplicación del método de Runge-Kutta a cada una de las ecuaciones diferenciales, con las condición inicial siguiente, en el instante t0

se esquematiza en la tabla adjunta. Como vemos además de los cuatro números k1, k2, k3, k4 para la primera ecuación diferencial precisamos otros cuatro números l1, l2, l3, l4 para la segunda ecuación diferencial. A partir del valor de x en el instante t, se determina el valor de x en el instante t+h, y a partir del valor de y en el instante t se determina el valor de y en el instante t+h mediante las fórmulas de la última fila de la tabla.

dx dt =f(t,x,y)

dy dt =g(t,x,y)

k 1 =h·f(t,x,y) k 2 =h·f( t+ 1 2 h,x+ 1 2 k 1 ,y+ 1 2 l 1 , ) k 3 =h·f( t+ 1 2 h,x+ 1 2 k 2 ,y+ 1 2 l 2 ) k 4 =h·f( t+h,x+ k 3 ,y+ l 3 )

l 1 =h·g(t,x,y) l 2 =h·g( t+ 1 2 h,x+ 1 2 k 1 ,y+ 1 2 l 1 ) l 3 =h·g( t+ 1 2 h,x+ 1 2 k 2 ,y+ 1 2 l 2 ) l 4 =h·g( t+h,x+ k 3 ,y+ l 3 )

x(t+h)=x(t)+ 1 6 ( k 1 +2 k 2 +2 k 3 + k 4 )

y(t+h)=y(t)+ 1 6 ( l 1 +2 l 2 +2 l 3 + l 4 )

Definimos la función rk_2_1 que resuelve el sistema de dos ecuaciones diferenciales de primer orden, cuando le pasamos:

Nos devuelve los vectores x e y para cada instante que se guarda en el vector t comprendido entre el instante inicial t0 y el final tf.

function [t,x,y] =rk_2_1(f,g,t0,tf,x0,y0,n)
      h=(tf-t0)/n;
      t=t0:h:tf;
%reserva memoria para n+1 element(i)os del vect(i)or x(i)
      x=zeros(n+1,1); 
      y=zeros(n+1,1);
      x(1)=x0; y(1)=y0;
      
      for i=1:n
        k1=h*f(t(i),x(i),y(i));
        l1=h*g(t(i),x(i),y(i));
        k2=h*f(t(i)+h/2,x(i)+k1/2,y(i)+l1/2);
        l2=h*g(t(i)+h/2,x(i)+k1/2,y(i)+l1/2);
        k3=h*f(t(i)+h/2,x(i)+k2/2,y(i)+l2/2);
        l3=h*g(t(i)+h/2,x(i)+k2/2,y(i)+l2/2);
        k4=h*f(t(i)+h,x(i)+k3,y(i)+l3);
        l4=h*g(t(i)+h,x(i)+k3,y(i)+l3);
        
        x(i+1)=x(i)+(k1+2*k2+2*k3+k4)/6;
        y(i+1)=y(i)+(l1+2*l2+2*l3+l4)/6;        
      end
end

Consideremos una serie radioactiva de tres elementos A-->B-->C en la que, una sustancia radiactiva A se desintegra y se transforma en otra sustancia radiactiva B, que a su vez se desintegra y se transforma en una sustancia C estable. Las ecuaciones diferenciales que gobiernan el proceso y sus soluciones analíticas son, respectivamente,

dx dt =axx= x 0 exp(at) dy dt =axbyy= a ba x 0 ( exp(at)exp(bt) )

La solución analítica que aparece a la derecha, se ha obtenido con las condiciones iniciales t=0, x=x0 e y=0. La segunda solución se obtiene siempre que a sea distinto de b. En el caso de que a sea igual a b, la solución analítica para y es

y= x 0 aexp(at)

La interpretación del sistema de ecuaciones diferenciales no es complicada. En la unidad de tiempo, desaparecen ax núcleos de la sustancia A al desintegrarse (primera ecuación). En la unidad de tiempo, se producen ax núcleos de la sustancia B y a su vez desaparecen bx núcleos de la sustancia B, que al desintegrarse se transforman en núcleos de la sustancia C estable (segunda ecuación).

Escribimos un script en el que definiremos las funciones f(t,x,y), g(t,x,y), las condiciones iniciales y llamaremos a la función rk_2_1

a=input('parámetro a: ');
b=input('parámetro b: ');
x0=input('valor inicial de x: ');
y0=input('valor inicial de y: ');
tf=input('tiempo final, tf: ');
n=input('número de pasos, n: ');
f=@(t,x,y) -a*x; 
g=@(t,x,y) a*x-b*y;
%condiciones iniciales
t0=0; 

[t,x,y]=rk_2_1(f,g,t0,tf,x0,y0,n);
hold on
plot(t,x,'b')
plot(t,y,'r')
xlabel('t')
ylabel('x,y')
grid on
legend('x(t)','y(t)')
title('dx/dt=-ax, dy/dt=ax-by')
hold off

En la ventana de comandos corremos el script

parámetro a: 0.1
parámetro b: .2
valor inicial de x: 100
valor inicial de y: 0
tiempo final, tf: 10
número de pasos, n: 40

Ecuación diferencial de segundo orden

Existen muchas situaciones en las que es necesario resolver una ecuación diferencial de segundo orden.

d 2 x d t 2 =f(t,x,v)

con las condiciones iniciales

x( t 0 )= x 0 ( dx dt ) t 0 = v 0

Una ecuación diferencial de segundo orden es equivalente a un sistema de dos ecuaciones diferenciales de primer orden, por lo que aplicaremos el mismo esquema.

dx dt =v

dv dt =f(t,x,v)

k 1 =hv k 2 =h( v+ 1 2 l 1 ) k 3 =h( v+ 1 2 l 2 ) k 4 =h( v+ l 3 )

l 1 =h·f(t,x,v) l 2 =h·f( t+ 1 2 h,x+ 1 2 k 1 ,v+ 1 2 l 1 ) l 3 =h·f( t+ 1 2 h,x+ 1 2 k 2 ,v+ 1 2 l 2 ) l 4 =h·f( t+h,x+ k 3 ,v+ l 3 )

x(t+h)=x(t)+ 1 6 ( k 1 +2 k 2 +2 k 3 + k 4 )

v(t+h)=v(t)+ 1 6 ( l 1 +2 l 2 +2 l 3 + l 4 )

Definimos la función rk_2 que resuelve la ecuación diferencial de segundo orden, cuando le pasamos:

Nos devuelve los vectores de las posiciones x y las velocidades v para cada instante que se guarda en el vector t comprendido entre el instante inicial t0 y el final tf.

function [t,x,v] =rk_2(f,t0,tf,x0,v0,n)
      h=(tf-t0)/n;
      t=t0:h:tf;
%reserva memoria para n+1 element(i)os del vect(i)or x(i)
      x=zeros(n+1,1); 
      v=zeros(n+1,1);
      x(1)=x0; v(1)=v0;
      
      for i=1:n
        k1=h*v(i);
        l1=h*f(t(i),x(i),v(i));
        k2=h*(v(i)+l1/2);
        l2=h*f(t(i)+h/2,x(i)+k1/2,v(i)+l1/2);
        k3=h*(v(i)+l2/2);
        l3=h*f(t(i)+h/2,x(i)+k2/2,v(i)+l2/2);
        k4=h*(v(i)+l3);
        l4=h*f(t(i)+h,x(i)+k3,v(i)+l3);
        
        x(i+1)=x(i)+(k1+2*k2+2*k3+k4)/6;
        v(i+1)=v(i)+(l1+2*l2+2*l3+l4)/6;        
      end
end

La ecuación diferencial que describe un oscilador armónico amortiguado y su solución para unas condiciones iniciales fijadas es

d 2 x d t 2 +2γ dx dt + ω 0 2 x=0 ω= ω 0 2 γ 2 x=exp(γt)( Asin(ωt)+Bcos(ωt) ) v= dx dt =γexp(γt)( Asin(ωt)+Bcos(ωt) )+ωexp(γt)( Acos(ωt)Bsin(ωt) )

Condiciones iniciales: en el instante t=0, la posición inicial es x0 y la velocidad inicial v0.

t=0{ x 0 =B v 0 =γB+ωA x=exp(γt)( v 0 +γ x 0 ω sin(ωt)+ x 0 cos(ωt) )

Escribimos un script en el que definiremos la función f(t,x,v), las condiciones iniciales y llamaremos a la función rk_2

w0=input('frecuencia angular w0: ');
g=input('rozamiento, gamma: ');
x0=input('posición inicial, x0: ');
v0=input('velocidad inicial,v0: ');
tf=input('tiempo final, tf: ');
n=input('número de pasos, n: ');
f=@(t,x,v) -2*g*v-w0*w0*x; 
%condiciones iniciales
t0=0;
hold on
%solución numérica
[t,x,v]=rk_2(f,t0,tf,x0,v0,n);
plot(t,x,'b')
%solución analítica
w=sqrt(w0*w0-g*g);
x=((v0+g*x0)*sin(w*t)/w+x0*cos(w*t)).*exp(-g*t);
plot(t,x,'r')
grid on
xlabel('t')
ylabel('x');
legend('aproximado','exacto')
title('oscilador amortiguado')
hold off

En la ventana de comandos corremos el script oscilador con distintas condiciones iniciales

>> oscilador
frecuencia angular, w0: 2
rozamiento, gamma: 0.5
posición inicial, x0: 1.5
velocidad inicial, v0: 0
tiempo final, tf: 8
número de pasos, n: 100

No se aprecia tampoco diferencia entre la solución exacta y la numérica, aplicando el procedimiento de Runge_Kutta.