# solve and plot a system of nonlinear 2nd order differential equations

16 views (last 30 days)
Zhen Zhen on 25 Mar 2019
Edited: Lewis Fer on 12 Jun 2021
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)
Cheers,
Jane
##### 2 CommentsShowHide 1 older comment
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
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!

Lewis Fer on 10 Jun 2021
Edited: Lewis Fer on 12 Jun 2021
Hello, I am having troubles solving a system of second order nonlinear equations with boundary conditions using MATALB
Here is the equations:
f''(t)=3*f(t)*g(t) +5;
g''(t)=-4f(t)*g(t)+7;
the boundary conditions are: f'(0)=0 et g'(o)=5;
g(0)=3 et f'(2)=f(2)