Asked by Ryan Bowman
on 6 Dec 2018 at 2:27

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

Answer by madhan ravi
on 6 Dec 2018 at 3:24

Edited by madhan ravi
on 6 Dec 2018 at 3:37

Accepted Answer

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
on 6 Dec 2018 at 3:51

madhan ravi
on 6 Dec 2018 at 4:40

Jesse Crotts
on 9 Dec 2018 at 0:12

I have a very similar problem at:

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.

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 4 Comments

## madhan ravi (view profile)

Direct link to this comment:https://www.mathworks.com/matlabcentral/answers/434127-how-to-solve-system-of-2nd-order-differential-equations-using-ode45#comment_646891

## Star Strider (view profile)

Direct link to this comment:https://www.mathworks.com/matlabcentral/answers/434127-how-to-solve-system-of-2nd-order-differential-equations-using-ode45#comment_646892

## Ryan Bowman (view profile)

Direct link to this comment:https://www.mathworks.com/matlabcentral/answers/434127-how-to-solve-system-of-2nd-order-differential-equations-using-ode45#comment_646893

## Ryan Bowman (view profile)

Direct link to this comment:https://www.mathworks.com/matlabcentral/answers/434127-how-to-solve-system-of-2nd-order-differential-equations-using-ode45#comment_646898

Sign in to comment.