How to use ODE45 for a coupled system of differential equations at a specific time point
12 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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1347369/image.png)
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1347574/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1347579/image.png)
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
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!