Issues solving a system of 2, second order, differential equations with ode45

Hello everybody,

I am starting with Matlab to solve a system of 2 differential equations, second order.

My code is below. It works, but the output is totally wrong.

This code is something taken on matlab help center and adapted to my problem, to help me start.

I have written things in the code "#comment " where it's not clear to me.

Your help would be very helpful to me.

Thanks a lot in advance.

clear all; 
close all;

syms x(t) y(t)

%% Constants  
Cd      = 0,45                                                             ;% drag coef aircraft          [-] 
Af      = 2                                                                ;% reference area aircraft     [-] 
m       = 10                                                               ;% total mass aircraft         [kg]
g       = 9.81                                                             ;% acc° of gravity             [m/s2]
Teta    = 85                                                               ;% launch angle                [deg]
V0      = 83                                                               ;% launch velocity             [m/s]
rho_Air = 1,2                                                              ;% air density                 [kg/m3]
t0      = 0;                                                               ;% time interval, initial t    [s]
tf      = 5400                                                             ;% time interval, final t      [s]


%% Initial Conditions  
x_to  = 0                                                                  ;% x(t=0)= 0                    [m] 
dx_t0 = V0*cosd(Teta)                                                      ;% x'(t=0)= V0*cosd(Teta)       [m/s] 
y_to  = 0                                                                  ;% y'(t=0)= 0                   [m] 
dy_t0 = V0*sind(Teta)                                                      ;% x(t=0)= 0                    [m/s]
y0    = [x_to dx_t0 y_to dy_t0 ]                                           ;% vector initial conditions 


%% System differential equations to be solved
eq1 = diff(x,2) == -( (Cd*Af*rho_Air)/(2*m) ) * diff(x,1) * (x^2 + y^2)^(1/2)
eq2 = diff(y,2) == - g - ( (Cd*Af*rho_Air)/(2*m) ) * diff(y,1) * ( (diff(x,1))^2 + (diff(y,1))^2 )^(1/2)

%% variables  
vars = [x(t); y(t)]

V = odeToVectorField([eq1,eq2])
M = matlabFunction(V,'vars', {'t','Y'}); 
% #comment: why the variable x doesn't appear here ? as I solve on x(t) and
% y(t) -> in V I call eq1 and eq2, then it's x(t) and y(t)

interval = [t0 tf];  %time interval    

ySol = ode45(M,interval,y0); 
% #comment: vector y0 containts the boundary conditions for both x and y. 
% As here I solve on y, why do we have boundary conditions in x as well ?
% is the function ode capable to make the difference in between boudady
% conditions for x and y ? 
tValues = linspace(interval(1),interval(2),1000);
yValues = deval(ySol,tValues,1); %number 1 denotes first position likewise you can mention 2 to 6 for the next 5 positions


plot(tValues,yValues)   %if you want to plot position vs time just swap here

% #comment: I need to solve on x as well
% thank you !  

Answer 1

[Torsten]

Replace

Cd      = 0,45                                                             ;% drag coef aircraft  [-]
rho_Air = 1,2                                                              ;% air density

by

 Cd      = 0.45                                                            ;% drag coef aircraft  [-] 
 rho_Air = 1.2                                                             ;% air density
 

If you use

 [V,S] = odeToVectorField([eq1,eq2])

instead of

 V = odeToVectorField([eq1,eq2]) 

you can see the order of the variables in the numerical solution vector Y comprising Y = [y,dy,x,dx].

You will notice that you have to change your initial conditions vector from

y0    = [x_to dx_t0 y_to dy_t0 ]                                           ;% vector initial conditions 

to

y0    = [ y_to dy_t0 x_to dx_t0]                                           ;% vector initial conditions 

[Labarrere]

Thank you Torsten.

It's better, but the final output is still wrong.

It feels like I don't use the good method to solves this system of equations.

Written properly the system of ode I need to solve looks like that:

Instead, I should use the proper ode45 solver.

% then equation 1 in the previous code becomes
f1 = @(t,x,u) -( (Cd*Af*rho_Air)/(2*m) ) * u * ( (x)^2 + (y)^2 )^(1/2);

% and boundary condiions are 
x(1)=0
u(1)=5


% then equation 2 in the previous code becomes
f2 = @(t,y,v) -g-( (Cd*Af*rho_Air)/(2*m) ) * v * ( (u)^2 + (v)^2 )^(1/2)
% and boundary condiions are 
t(1)=0
y(1)=0
v(1)=10

If you have any tips, please don't hesitate, will be much appreciated

[Torsten]

There was an error in eq1. Now the equations are set according to your physical model. If the results (I plotted y) are still wrong, you will have to check your equations and/or your parameters.

clear all; 
close all;
syms x(t) y(t)

%% Constants  
Cd      = 0.45                                                             ;% drag coef aircraft          [-] 
Af      = 2                                                                ;% reference area aircraft     [-] 
m       = 10                                                               ;% total mass aircraft         [kg]
g       = 9.81                                                             ;% acc° of gravity             [m/s2]
rho_Air = 1.2                                                              ;% air density                 [kg/m3]
t0      = 0;                                                               ;% time interval, initial t    [s]
tf      = 5400                                                             ;% time interval, final t      [s]
%% Initial Conditions  
x_to  = 0                                                                  ;% x(t=0)= 0                    [m] 
dx_t0 = 5                                                                  ;% x'(t=0)= V0*cosd(Teta)       [m/s] 
y_to  = 0                                                                  ;% y'(t=0)= 0                   [m] 
dy_t0 = 10                                                                 ;% x(t=0)= 0                    [m/s]
y0    = [y_to dy_t0 x_to dx_t0  ]                                          ;% vector initial conditions 
%% System differential equations to be solved
eq1 = diff(x,2) == - Cd*Af*rho_Air/(2*m)*diff(x)*(diff(x)^2 + diff(y)^2)^(1/2)

eq1(t) = 

eq2 = diff(y,2) == - g - Cd*Af*rho_Air/(2*m)*diff(y)*(diff(x)^2 + diff(y)^2)^(1/2)

eq2(t) = 

%% variables  
%vars = [x(t); y(t)]
[V,S] = odeToVectorField([eq1,eq2])

V = 

S = 

●

M = matlabFunction(V,'vars', {'t','Y'}); 
% #comment: why the variable x doesn't appear here ? as I solve on x(t) and
% y(t) -> in V I call eq1 and eq2, then it's x(t) and y(t)
interval = [t0 tf];  %time interval    
ySol = ode45(M,interval,y0); 
tValues = linspace(interval(1),interval(2),1000);
yValues = deval(ySol,tValues,1); %number 1 denotes first position likewise you can mention 2 to 6 for the next 5 positions
plot(tValues,yValues)   %if you want to plot position vs time just swap here

[Torsten]

If you use "ode45" directly, the code would be

clear all 
close all
%% Constants  
Cd      = 0.45                                                             ;% drag coef aircraft          [-] 
Af      = 2                                                                ;% reference area aircraft     [-] 
m       = 10                                                               ;% total mass aircraft         [kg]
g       = 9.81                                                             ;% acc° of gravity             [m/s2]
rho_Air = 1.2                                                              ;% air density                 [kg/m3]
t0      = 0                                                                ;% time interval, initial t    [s]
tf      = 5400                                                             ;% time interval, final t      [s]
%% Initial Conditions  
x_to  = 0                                                                  ;% x(t=0)= 0                    [m] 
dx_t0 = 5                                                                  ;% x'(t=0)= V0*cosd(Teta)       [m/s] 
y_to  = 0                                                                  ;% y'(t=0)= 0                   [m] 
dy_t0 = 10                                                                 ;% x(t=0)= 0                    [m/s]
y0    = [ x_to dx_t0 y_to dy_t0 ]                                          ;% vector initial conditions 
% Equations
fun = @(t,y)[y(2);...
             - Cd*Af*rho_Air/(2*m)*y(2)*(y(2)^2 + y(4)^2)^(1/2);...
             y(4);
             -g - Cd*Af*rho_Air/(2*m)*y(4)*(y(2)^2 + y(4)^2)^(1/2)];
% Integration interval
tspan = [t0 tf];  %time interval
% Call solver
[T,Y] = ode45(fun,tspan,y0);
% Plot solutions
plot(T,Y(:,3))

●