Uso de Runge-Kutta de cuarto orden para resolver una oda de segundo orden en MATLAB

Question

Para una bola de masa $m$ , radio $r$ cayendo libremente a través del aire de densidad $p$ La fuerza de resistencia del aire es proporcional a la densidad del aire y al cuadrado de la velocidad de la pelota, $v$ y el área de una proyección de la bola sobre un plano perpendicular a la dirección de su movimiento. Así, en una dimensión, la ecuación del movimiento es:

$mx'' = mg - 0.5CApv^2$ donde g es la aceleración debida a la gravedad

Si en cambio se considera el movimiento en dos dimensiones, entonces las siguientes ecuaciones lo describen:

$mx'' = -0.5CApv_xv$

$my'' = mg - 0.5CApv_yv$

Con v \= $(v_x,v_y)$ y $v = (v_x^2 + v_y^2)$

Ahora quiero resolver las dos ecuaciones diferenciales de segundo orden anteriores utilizando el método Runge-Kutta de cuarto orden en MATLAB. Utilicé el siguiente bolo de código:

SOLUCIÓN EN MATLAB

%NOTES:

%Fouth order Runge-Kutta to solve two systems of second order diff
%equations:
%  x'' = -(1/2CAp/m)v_xV = -Bv_xV ; y'' = g - Bv_yV

%For x: use substitutions: x_1 =x ; x_2 = x_1' = v_x
%for y: use : y_1 = y ; y_2 = y_1' = v_y

%Now we have these two pairs of 1st order equations:

%1) x_1' = x_2 = f(x_1,x_2,y_1,y_2);
%2) x_2' = -Bx_2sqrt(x_2^2 + y_2^2) = g(x_1,x_2,y_1,y_2);
%3) y_1' = y_2 = h(x_1,x_2,y_1,y_2);    
%4) y_2' = g - By_2sqrt(x_2^2 + y_2^2) = q(x_1,x_2,y_1,y_2);

g = 9.81;   %acceleration due to gravity
C = 0.5;    %a constant C, set to 0.5
p = 1.225;  %the density of air (in SI units)
m = 2.47e-03;  %the mass of the ball (SI)
r = 20e-03;    %the radius of the ball (SI)
A = pi*(r^2);  %the projection of the ball onto a normal plane, which is a circle of radius r
B = (C*A*p)/(2*m);  %a combination of all the constants
t_stop = 10;    %the end time(how long I want to integrate the motion)
N = 500;        %number of iterations
dt = t_stop/N;     %time step

x_1(1) = 0;         %initializing the x_1 - arrray
x_2(1) = 7*cos(pi/4);   %initializing the x_2 - arrray
y_1(1) = 0;     %initializing the y_1 - arrray
y_2(1) = 7*sin(pi/4);   %initializing the y_2 - arrray
t(1) = 0;       %the initial time

for n = 1:N;    
    t(n+1) = n*dt;

    k1 = dt*x_2(n);
    m1 = dt*-B*x_2(n)*sqrt(x_2(n)^2 + y_2(n)^2);
    L1 = dt*(y_2(n));
    Z1 = dt*(9.81 - B*y_2(n)*sqrt(x_2(n)^2 + y_2(n)^2));

    k2 = dt*(x_2(n)+m1/2);
    m2 = dt*-B*(x_2(n)+ m1/2)*sqrt((x_2(n)+ m1/2)^2 + (y_2(n)+Z1/2)^2);
    L2 = dt*(y_2(n)+ Z1/2);
    Z2 = dt*(9.81 - B*(y_2(n)+Z1/2)*sqrt((x_2(n)+m1/2)^2 + (y_2(n)+Z1/2)^2));

    k3 = dt*(x_2(n)+m2/2);
    m3 = dt* -B*(x_2(n)+ m2/2)*sqrt((x_2(n)+ m2/2)^2 + (y_2(n)+Z2/2)^2);
    L3 = dt*(y_2(n)+ Z2/2);
    Z3 = dt*(9.81 - B*(y_2(n)+Z2/2)*sqrt((x_2(n)+m2/2)^2 + (y_2(n)+Z2/2)^2));    

    k4 = dt*(x_2(n)+ m3);
    m4 = dt*-B*(x_2(n)+ m3)*sqrt((x_2(n)+ m3)^2 + (y_2(n)+Z3)^2);
    L4 = dt*(y_2(n)+ Z3);
    Z4 = dt*(9.81 - B*(y_2(n)+Z3)*sqrt((x_2(n)+ m3)^2 + (y_2(n)+Z3)^2));   

    x_1(n+1) = x_1(n)+1/6*(k1 + 2*k2 + 2*k3 + k4);
    x_2(n+1) = x_2(n)+1/6*(m1 + 2*m2 + 2*m3 + m4);
    y_1(n+1) = y_1(n)+1/6*(L1 + 2*L2 + 2*L3 + L4);
    y_2(n+1) = y_2(n)+1/6*(Z1 + 2*Z2 + 2*Z3 + Z4);

end

plot(y_1,x_1)

Problema : Cuando ejecuto el programa anterior, me da un gráfico bastante extraño que no se corresponde en absoluto con el movimiento de un proyectil. ¿Alguien puede ver cuál es el problema en mi código?

Answer 1

Tus gráficos se ven bien excepto que has definido y como positivo hacia abajo usando y"=g + etc. Si quieres que la y sea positiva hacia arriba, tienes que usar y"=g + etc.

El código ordenado con los gráficos de los componentes individuales del estado contra el tiempo, y el positivo de y hacia arriba es:

    %NOTES:

%Fouth order Runge-Kutta to solve two systems of second order diff
%equations:
%  x'' = -(1/2CAp/m)v_xV = -Bv_xV ; y'' = g - Bv_yV

%For x: use substitutions: x_1 =x ; x_2 = x_1' = v_x
%for y: use : y_1 = y ; y_2 = y_1' = v_y

%Now we have these two pairs of 1st order equations:

%1) x_1' = x_2 = f(x_1,x_2,y_1,y_2);
%2) x_2' = -Bx_2sqrt(x_2^2 + y_2^2) = g(x_1,x_2,y_1,y_2);
%3) y_1' = y_2 = h(x_1,x_2,y_1,y_2);    
%4) y_2' = g - By_2sqrt(x_2^2 + y_2^2) = q(x_1,x_2,y_1,y_2);

g = -9.81;   %acceleration due to gravity
C = 0.5;    %a constant C, set to 0.5
p = 1.225;  %the density of air (in SI units)
m = 2.47e-03;  %the mass of the ball (SI)
r = 20e-03;    %the radius of the ball (SI)
A = pi*(r^2);  %the projection of the ball onto a normal plane, which is a circle of radius r
B = (C*A*p)/(2*m);  %a combination of all the constants
t_stop = 10;    %the end time(how long I want to integrate the motion)
N = 500;        %number of iterations
dt = t_stop/N;     %time step

x_1(1) = 0;         %initializing the x_1 - arrray
x_2(1) = 7*cos(pi/4);   %initializing the x_2 - arrray
y_1(1) = 0;     %initializing the y_1 - arrray
y_2(1) = 7*sin(pi/4);   %initializing the y_2 - arrray
t(1) = 0;       %the initial time

for n = 1:N;    
    t(n+1) = n*dt;

    k1 = dt*x_2(n);
    m1 = dt*-B*x_2(n)*sqrt(x_2(n)^2 + y_2(n)^2);
    L1 = dt*(y_2(n));
    Z1 = dt*(g - B*y_2(n)*sqrt(x_2(n)^2 + y_2(n)^2));

    k2 = dt*(x_2(n)+m1/2);
    m2 = dt*-B*(x_2(n)+ m1/2)*sqrt((x_2(n)+ m1/2)^2 + (y_2(n)+Z1/2)^2);
    L2 = dt*(y_2(n)+ Z1/2);
    Z2 = dt*(g - B*(y_2(n)+Z1/2)*sqrt((x_2(n)+m1/2)^2 + (y_2(n)+Z1/2)^2));

    k3 = dt*(x_2(n)+m2/2);
    m3 = dt* -B*(x_2(n)+ m2/2)*sqrt((x_2(n)+ m2/2)^2 + (y_2(n)+Z2/2)^2);
    L3 = dt*(y_2(n)+ Z2/2);
    Z3 = dt*(g - B*(y_2(n)+Z2/2)*sqrt((x_2(n)+m2/2)^2 + (y_2(n)+Z2/2)^2));    

    k4 = dt*(x_2(n)+ m3);
    m4 = dt*-B*(x_2(n)+ m3)*sqrt((x_2(n)+ m3)^2 + (y_2(n)+Z3)^2);
    L4 = dt*(y_2(n)+ Z3);
    Z4 = dt*(g - B*(y_2(n)+Z3)*sqrt((x_2(n)+ m3)^2 + (y_2(n)+Z3)^2));   

    x_1(n+1) = x_1(n)+1/6*(k1 + 2*k2 + 2*k3 + k4);
    x_2(n+1) = x_2(n)+1/6*(m1 + 2*m2 + 2*m3 + m4);
    y_1(n+1) = y_1(n)+1/6*(L1 + 2*L2 + 2*L3 + L4);
    y_2(n+1) = y_2(n)+1/6*(Z1 + 2*Z2 + 2*Z3 + Z4);

end

subplot(2,2,1);plot(t,x_1);title('x against time');
subplot(2,2,2);plot(t,x_2);title('v_x against time');

subplot(2,2,3);plot(t,y_1);title('y against time');
subplot(2,2,4);plot(t,y_2);title('v_y against time');

Que produce parcelas:

a partir de la cual vemos que v_x cae a cero, y v_y se aproxima a una velocidad terminal.