How to use ODE45 for a coupled system of differential equations at a specific time point
3 views (last 30 days)
Show older comments
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
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.
Accepted Answer
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)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!