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

12 Ansichten (letzte 30 Tage)
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 Kommentare
Sam Chak
Sam Chak am 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.

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Jon
Jon am 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

Weitere Antworten (0)

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Produkte


Version

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by