Clear Filters
Clear Filters

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

3 views (last 30 days)
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.
  3 Comments
Sam Chak
Sam Chak on 6 Apr 2023
@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.

Sign in to comment.

Accepted Answer

Jon
Jon on 6 Apr 2023
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

More Answers (0)

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!