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

147 views (last 30 days)

Show older comments

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

### Accepted Answer

madhan ravi
on 6 Dec 2018

Edited: madhan ravi
on 6 Dec 2018

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:

##### 5 Comments

Shantanu Chhaparia
on 20 Feb 2022

### More Answers (1)

Bob
on 14 Feb 2023

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

##### 0 Comments

### See Also

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!