MATLAB: Solve and plot a system of nonlinear 2nd order differential equations

differential equationssystem of differential equations

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

Best 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
Related Question