Asked by Zhen Zhen
on 25 Mar 2019

hi there,

I'm trying to plot a graph of against with the following equations of motion:

I've tried dsolve and ode45 yet there always seems to be some problems. I think ode45 might work better because apparently it would be easier to plot the graph by using some numerical method?

Here's my failed attempt to solve it: (I just set some variables to be 1 to make the problem easier here)

syms theta(t) phi(t) psi(t) C

dtheta = diff(theta , t);

d2theta = diff(theta , t , 2);

dphi = diff(phi , t);

%dpsi = diff(psi , t);

%alpha = C * (dpsi + dphi * cos(theta));

%beta = alpha * cos(theta) + dphi * (sin(theta))^2;

alpha = 1;

beta = 1;

dpsi = 1;

eqn1 = dphi == (beta - alpha * cos(theta)) / (sin(theta))^2 ; %equations of motion

eqn2 = d2theta == (dphi*(dphi * cos(theta) - alpha)+1)*sin(theta) ; %equations of motion

eqns = [eqn1 , eqn2];

%cond = [dpsi == 1];

[thetaSol(t) phiSol(t)] = dsolve(eqns)

Thanks a lot for your help and time in advance!

Cheers,

Jane

Answer by Teja Muppirala
on 26 Mar 2019

Edited by Teja Muppirala
on 26 Mar 2019

Accepted Answer

If you just need a plot and not a closed-form solution, then I'd recommend just using ODE45 without worrying about symbolic stuff. This is an example of how to solve this using ODE45 for initial conditions psi(0) = 0, theta(0) = 0, thetadot(0) = 1 over the time span [0 10].

function doODE

a = 1; % alpha

b = 1; % beta

% Define variables as:

% Y(1): phi

% Y(2): theta

% Y(3): thetadot

ic = [0;0;1]; % Pick some initial Conditions

tspan = [0 10];

[tout, yout] = ode45(@deriv,tspan,ic);

subplot(1,3,1), plot(tout,yout(:,1)); title('psi(t)')

subplot(1,3,2), plot(tout,yout(:,2)); title('theta(t)')

subplot(1,3,3), plot(yout(:,1),yout(:,2)); title('theta vs psi')

function dY = deriv(t,Y)

dY = [dpsi(Y(2)); ...

Y(3); ...

[dpsi(Y(2))*(dpsi(Y(2))*cos(Y(2))-a)+1]*sin(Y(2)) ];

end

function dp = dpsi(th)

dp = (b - a*cos(th))./(sin(th)).^2 ;

if isnan(dp) % Need to take care at theta = 0

dp = 0.5;

end

end

end

Zhen Zhen
on 28 Mar 2019

Thank you very much! it worked well!

I was wondering, if now I need to use the real alpha and beta (alpha = C * (dpsi + dphi * cos(theta)); beta = alpha * cos(theta) + dphi * (sin(theta))^2;) I know they are constants. Where should I put them in the code?

And if I need to print one more constant (say, E = 0.5 * (Y(3))^2 + (Y(1))^2 * (sin(Y(2)))^2) +0.5* a^2 + cos(Y(2)) ) where should I put it?

Thanks so much for your help!

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 2 Comments

## Walter Roberson (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/452333-solve-and-plot-a-system-of-nonlinear-2nd-order-differential-equations#comment_685783

## Walter Roberson (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/452333-solve-and-plot-a-system-of-nonlinear-2nd-order-differential-equations#comment_685784

Sign in to comment.