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

Sistema de dos ecuaciones diferenciales de segundo orden

El procedimiento numérico de Runge-Kutta puede aplicarse para resolver un sistema de dos ecuaciones diferenciales de segundo orden.

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

con las condiciones iniciales

x( t 0 )= x 0 ( dx dt ) t 0 = v x0 y( t 0 )= y 0 ( dy dt ) t 0 = v y0

Este sistema, se puede transformar en un sistema equivalente formado por cuatro ecuaciones diferenciales de primer orden. Aplicando dos veces el esquema descrito para una ecuación diferencial de segundo orden, obtenemos el esquema descrito en la siguiente tabla

dx dt = v x

d v x dt =f(t,x, v x ,y, v y )

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

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

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

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

dy dt = v y

d v y dt =g(t,x, v x ,y, v y )

q 1 =h v y q 2 =h( v y + 1 2 m 1 ) q 3 =h( v y + 1 2 m 2 ) q 4 =h( v y + m 3 )

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

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

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

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

Nos devuelve los vectores posición (x,y) y velocidad (vx,vy) para cada instante que se guarda en el vector t comprendido entre el instante inicial t0 y el final tf.

function [t,x,vx,y,vy]=rk_2_2(f,g,t0,tf,x0,vx0,y0,vy0,n)
      h=(tf-t0)/n;
      t=t0:h:tf;
      x=zeros(n+1,1); 
%reserva memoria para n+1 element(i)os del vect(i)or x(i)
      y=zeros(n+1,1);
      vx=zeros(n+1,1);
      vy=zeros(n+1,1);
      x(1)=x0; vx(1)=vx0;
      y(1)=y0; vy(1)=vy0;
  
     for i=1:n
        k1=h*vx(i);
        l1=h*f(t(i),x(i),vx(i),y(i),vy(i));
        q1=h*vy(i);
        m1=h*g(t(i),x(i),vx(i),y(i),vy(i));

        k2=h*(vx(i)+l1/2);
        l2=h*f(t(i)+h/2,x(i)+k1/2,vx(i)+l1/2,y(i)+q1/2,vy(i)+m1/2);
        q2=h*(vy(i)+m1/2);
        m2=h*g(t(i)+h/2,x(i)+k1/2,vx(i)+l1/2,y(i)+q1/2,vy(i)+m1/2);

        k3=h*(vx(i)+l2/2);
        l3=h*f(t(i)+h/2,x(i)+k2/2,vx(i)+l2/2,y(i)+q2/2,vy(i)+m2/2);
        q3=h*(vy(i)+m2/2);
        m3=h*g(t(i)+h/2,x(i)+k2/2,vx(i)+l2/2,y(i)+q2/2,vy(i)+m2/2);

        k4=h*(vx(i)+l3);
        l4=h*f(t(i)+h,x(i)+k3,vx(i)+l3,y(i)+q3,vy(i)+m3);
        q4=h*(vy(i)+m3);
        m4=h*g(t(i)+h,x(i)+k3,vx(i)+l3,y(i)+q3,vy(i)+m3);

        x(i+1)=x(i)+(k1+2*k2+2*k3+k4)/6;
        vx(i+1)=vx(i)+(l1+2*l2+2*l3+l4)/6;
        y(i+1)=y(i)+(q1+2*q2+2*q3+q4)/6;
        vy(i+1)=vy(i)+(m1+2*m2+2*m3+m4)/6;        
    end
end

Uno de los ejemplos más interesantes de resolución de un sistema de ecuaciones diferenciales de segundo orden es la descripción del movimiento de los cuerpos celestes, el cual tiene una solución analítica sencilla en coordenadas polares. La trayectoria seguida por un planeta es una cónica, una elipse en particular, en uno de cuyos focos está el centro fijo de fuerzas, el Sol. En la figura, se muestra la fuerza que ejerce el Sol sobre un planeta, inversamente proporcional al cuadrado de las distancias que separan sus centros y también, se muestran sus componentes rectangulares

Teniendo en cuenta que la fuerza que ejerce el Sol sobre un planeta viene descrita por la ley de la Gravitación Universal

F=G Mm r 2

donde M es la masa del Sol, m la masa del planeta y r la distancia entre el centro del Sol y del planeta. Las componentes de la aceleración del planeta serán

a x = F x m = Fcosθ m =G M r 3 x a y = F y m = Fsinθ m =G M r 3 y

Uno de los problemas del tratamiento numérico con ordenador, es la de reducir el problema a números simples e inteligibles por el usuario de un vistazo. Las masa de los planetas y del Sol son números muy grandes: la masa de la Tierra es 5.98·1024 kg, y 1.98·1030 kg, la del Sol. La distancia media entre la Tierra y el Sol es también muy grande 1.49·1011 m, y la constante G es muy pequeña 6.67·10-11 en el Sistema Internacional de Unidades.

Escalas

Antes de resolver el sistema de ecuaciones diferenciales por procedimientos numéricos, es conveniente prepararlas para que el ordenador no maneje números excesivamente grandes o pequeños.

Establecemos un sistema de unidades en el que la longitud se mide en unidades astronómicas, la distancia media entre el Sol y la Tierra. L=una UA=1.496·1011 m y el tiempo en unidades de año, P=un año= 365.26 días=3.156·107 s.

En el nuevo sistema de unidades x=Lx’, t=P·t’, la primera ecuación diferencial se escribe

d 2 x ' dt ' 2 L P 2 = GM (x ' 2 +y ' 2 ) 3/2 x' L L 3 d 2 x ' dt ' 2 = GM P 2 L 3 x' (x ' 2 +y ' 2 ) 3/2

Como L es el semieje mayor de la órbita de la Tierra alrededor del Sol, P es el periodo o tiempo que tarda en dar una vuelta completa, y M es la masa del Sol. Por la tercera ley de Kepler, el término

GM P 2 L 3 =4 π 2

Volviendo a la notación x e y para la posición y t para el tiempo en el nuevo sistema de unidades. El sistema de ecuaciones diferenciales se escribe

d 2 x d t 2 =4 π 2 x ( x 2 + y 2 ) 3/2 d 2 y d t 2 =4 π 2 y ( x 2 + y 2 ) 3/2

Se resuelve por procedimientos numéricos con las condiciones iniciales siguientes: en el instante t=0, x=x0, y=0, vx=0, vy=v0y.

Principio de conservación de la energía

La energía total de la partícula es una constante del movimiento. La energía de la partícula de masa m en el instante inicial t=0 es

E 0 = 1 2 m v 0 2 GmM r 0

Cuando E0<0 la partícula permanece confinada en el espacio que rodea a los dos cuerpos. Cuando E0≥0 la partícula escapa al infinito

La energía de la partícula en el instante t es igual a

E= 1 2 m( v x 2 + v y 2 ) GmM x 2 + y 2

En el nuevo sistema de unidades establecido

v=v’·L/P, x=x’·L, y=y’·L, d=d’·L

E m = 1 2 v ' 2 L 2 P 2 G M 1 L 1 x ' 2 +y ' 2 E m = 1 2 v ' 2 L 2 P 2 4 π 2 L 2 P 2 1 x ' 2 +y ' 2 e ' = E m P 2 L 2 = 1 2 v ' 2 4 π 2 1 x ' 2 +y ' 2

Volviendo a la notación previa. Definimos una nueva energía e por unidad de masa en este sistema de unidades

e= 1 2 ( v x 2 + v y 2 )4 π 2 1 x 2 + y 2

El programa evalúa el cociente

| e e 0 e 0 |·100

que denominaremos tanto por ciento de error. Cuando la energía e difiere apreciablemente de e0, la trayectoria calculada puede que se desvíe significativamente de la real, y esto suele ocurrir cuando la partícula está muy cerca del centro de fuerzas

Elaboramos un script para describir el movimiento de un cuerpo celeste.

x0=input('posición inicial x: ');
vx0=input('velocidad incial x: ');
y0=0;
vy0=input('velocidad incial y: ');
tf=input('tiempo final, tf: ');
n=input('número de pasos, n: ');

f=@(t,x,vx,y,vy) -4*pi*pi*x/(sqrt(x*x+y*y))^3; 
g=@(t,x,vx,y,vy) -4*pi*pi*y/(sqrt(x*x+y*y))^3;
%condiciones iniciales
t0=0; 
e0=(vx0*vx0+vy0*vy0)/2-4*pi*pi/sqrt(x0*x0+y0*y0); %energía inicial

[t,x,vx,y,vy]=rk_2_2(f,g,t0,tf,x0,vx0,y0,vy0,n);
plot(x,y,'r')
xlabel('x')
ylabel('y');
grid on
axis equal
title('órbita de un planeta')
%energía final
ef=(vx(end)*vx(end)+vy(end)*vy(end))/2-4*pi*pi/sqrt(x(end)*x(end)+y(end)*y(end));
%error relativo
fprintf('Error en energía %1.4f\n',abs((ef-e0)/e0))

En la ventana de comandos corremos el script

posición inicial x: 1
velocidad incial x: 0
velocidad incial y: 6.27
tiempo final, tf: 1
número de pasos, n: 40
Error en energía 0.0000