How to use ODE45 for a coupled system of differential equations at a specific time point

Hey everyone. I have been using ode45 for a while and it was smooth for me. However, I am having trouble doing it for this system because its a bit complicated. I had multiple attempts fail before giving up. I hope someone can help.

I have the system shown below in the picture, its a system of three differential equations. Wx, Wy, and Wz are constants that I define, and the initial conditions for theta, phi, and psi, are also constants that I define. I want to find the solution for theta, phi and psi at a specific point of time only. Example, time = 2.2 seconds. Note that the dots above the variables indicate the derivative, so x_dot is dx/dt.

[Davide Masiello]

If you want to know the solution at t = 2.2 s you ought to integrate the system up to that time, there's no way to otherwise know it without integrating the ODEs.

[Sam Chak]

@Ali Almakhmari, Can you show your attempted code? We will investigate the problem why you cannot get the solution exactly at sec. However, a quick analysis on this set of ODEs shows that singularity can occur at rad.

[Ali Almakhmari]

Thank you so much guys. Your notes helped me a lot!.

Answer 1

[Jon]

So, expanding on @Davide Masiello comment, if all you want is the final solution at t = 2.2s you could wrap your call to ode45 using something like this (numerical values are just for illustrative purposes)

theta0 = pi/10;
phi0= pi/3;
psi0 = pi/9;
w = randn(3,1);
tFinal = 2.2
[theta,phi,psi] = solveSys(tFinal,theta0,phi0,psi0,w)

where solveSys is a function that you define (store as its own .m file )

function [theta,phi,psi] = solveSys(tFinal,theta0,phi0,psi0,w)
% solve system of odes to get final value of angles theta,phi and psi for
% specified w = [wx,wy,wz]
tspan = [0,tFinal];
x0 = [theta0,phi0,psi0]; % initial conditions
        % define some local variables for readability
        wx = w(1);
        wy = w(2);
        wz = w(3);
[~,x] = ode45(@fun,tspan,x0); % just get back angle values, not interested in times
% just return the final values
theta = x(end,1); 
phi = x(end,2);
psi = x(end,3);
    function xdot = fun(~,x)
        % calculates state derivatives, time is not needed, but needs to be
        % included in signature, so use ~
        % use nested function so parameters wx,wy,wz will be
        % available

        % define states for clearer readability
        theta = x(1);
        phi = x(2);
        psi = x(3);
        phiDot= wx + tan(theta)*sin(phi)*wy + tan(theta)*cos(phi)*wz;
        thetaDot = cos(phi)*wy - sin(phi)*wz;
        psiDot = sec(theta)*sin(phi)*wy + sec(theta)*cos(phi)*wz;

        xdot = [phiDot;thetaDot;psiDot];
    end
end

[Jon]

You're welcome