How to solve system of 2nd order differential equations using ode45

I have three 2nd order differential equations with my initial conditions and I'm trying to use the ode45 function in matlab to solve this. I wish to get the solution where my output is x,y,z position vs. time plot(2nd derivative) as well as a dx,dy,dz velocity vs. time plot. I get multiple errors and I'm not sure how to fix it. Here is my code:

%Clohessy-Wiltshire Equations
% d2x = 2*n*dy + 3*(n^2)*x;
% d2y = -2*n*dx;
% d2z = (-n^2)*z;
% 
% %Initial Conditions
% x(0) = -0.066538073651029; %km
% y(0) = 0.186268907590665; %km
% z(0) = 0.000003725378152; %km
% dx(0) = -0.000052436200437; %km/s
% dy(0) = 0.000154811363681; %km/s
% dz(0) = 0.000210975508077; %km/s

%Constants
a = 6793.137; %km
mu = 398600.5; %km^3/s^2
n = sqrt(mu/a^3);

t0 = 0; %seconds
tf = 5400; %seconds
initial_condition = [x(0) y(0) z(0)]

[T,position] = ode45(@(t,position)Clohessy_Wiltshire(t,x,y,z),[t0 tf],initial_condition);
figure(1), hold on, plot(T,position(:,1),'b'), plot(T,position(:,2), 'r'), plot(T,position(:,3), 'g')
title('Position(X,Y,Z) vs. Time')
ylabel('Position(X,Y,Z)(km)')
xlabel('Time(s)')

figure(2), hold on, plot(T,position(:,4),'b'), plot(T,position(:,5), 'r'), plot(T,position(:,6), 'g')
title('Velocity(dX,dY,dZ) vs. Time')
ylabel('Velocity(dX,dY,dZ)')
xlabel('Time(s)')

function position = Clohessy_Wiltshire(~,x,y,z,dx,dy,dz,n)
x(0) = -0.066538073651029;
dx(0) = -0.000052436200437;
dx(2) = 2*n*dy;
y(0) = 0.186268907590665;
dy(0) =  0.000154811363681;
dy(2) = -2*n*dx;
z(0) = 0.000003725378152;
dz(0) = 0.000210975508077;
dz(2) = -n^2*z;
end

[madhan ravi]

please upload the equation in latex form

[Star Strider]

Note that your ‘Clohessy_Wiltshire’ anonymous function never creates a variable named ‘position’.

It should return a column vector of derivatives.

[Ryan Bowman]

I'm not sure how to do that but these are the equations/initial conditions i'm given. Thanks.

[Ryan Bowman]

I'm essentially just integrating these 3 2nd order differential equations.

Answer 1

[madhan ravi]

Here you go!

syms x(t) y(t) z(t)
%Clohessy-Wiltshire Equations
% d2x = 2*n*dy + 3*(n^2)*x;
% d2y = -2*n*dx;
% d2z = (-n^2)*z;
%Constants
a = 6793.137; %km
mu = 398600.5; %km^3/s^2
n = sqrt(mu/a^3);
t0 = 0; %seconds
tf = 5400; %seconds
dx=diff(x,t);
dy=diff(y,t);
dz=diff(z,t);
%Initial Conditions
c1 = -0.066538073651029; %km
c2 =0.186268907590665; %km
c3 =0.000003725378152; %km
c4 = -0.000052436200437; %km/s
c5 =0.000154811363681; %km/s
c6 = 0.000210975508077; %km/s
y0 = [c1 c2 c3 c4 c5 c6];
eq1 = diff(x,2) == 2*n*dy + 3*(n^2)*x;
eq2 = diff(y,2) == -2*n*dx;
eq3 = diff(z,2) == (-n^2)*z;
vars = [x(t); y(t); z(t)]
V = odeToVectorField([eq1,eq2,eq3])
M = matlabFunction(V,'vars', {'t','Y'});
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

The graph of the first position looks like below:

[Ryan Bowman]

Thank you so much! Am I also able to get the x and z values on the same plot with this? I should have three different lines on this plot.

[madhan ravi]

Anytime :) , yes you can just use hold on and then plot the rest , if my answer helped make sure to accept the answer and give a vote.

[Jesse Crotts]

I have a very similar problem at:

https://www.mathworks.com/matlabcentral/answers/434640-how-to-simulate-using-ode23-with-multidimentional-second-order-differential-equations

I'm basically trying to use ode23 or ode45 to solve a system of 2nd order differential equations that look like this:

[M]*xdotdot+[K]*x=[Q]

where M and K are 10x10 matrices and Q is a 10x1 matrix.

I would appreciate any help with it.

[James VanderVeen]

You specifiy the variable:

vars = [x(t); y(t); z(t)]

But you never use it. Is there a reason for inserting this into the code?

[Shantanu Chhaparia]

hey! are there some examples (of system of higher order differential equation) on matlab site? If yes, can you please share the link. I was unable to locate them and I had certain doubts which might get cleared by looking over those. Thanks.

Answer 2

[Bob]

Hopefully, it is valid

% define n, where Earth GM : μ = 398600.442 km³/s²
n = sqrt(398600.442e9/earthRadius^3) ;  % note: earthRadius < a
% and matrices A and B
A = [0 2 0; -2 0 0; 0 0  0].* n ;
B = [3 0 0;  0 0 0; 0 0 -1].* n^2 ;
% Define symbolic variable t and vector u(t) ≡ [x; ẋ]
syms t ;
syms u(t) [6,1] matrix ;
% Specify initial value / start position
s = [-66.538073651029;186.268907590665;0.003725378152;...   % m
      -0.052436200437;  0.154811363681;0.210975508077]  ;   % m/s
% Define ODE function...
M = @(t,u)[u(4:6); A*u(4:6) + B*u(1:3)] ;
% ...and solve by ode45 on (0;5400] time interval, s = u(0)
z = ode45(M,[0 5400],s) ;
r = 0:54:5400 ;     % points range to plot the results
plot(r,deval(z,r,[1:3])); % distance [m] vs time

plot(r,deval(z,r,[4:6])); % velocity [m/s] vs time