48 views (last 30 days)

Show older comments

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

Walter Roberson
on 26 Mar 2019

Have a look at odeFunction() and in particular the first example. It shows you the functions you need to use in order to convert the symbolic forms into something you can call with ode45.

Use the options structure created by odeset to designate an OutputFcn . You might want to use @odephas2 to construct a 2D phase plot, perhaps.

Teja Muppirala
on 26 Mar 2019

Edited: Teja Muppirala
on 26 Mar 2019

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

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

Start Hunting!