# Time-based integration of ODEs

10 vues (au cours des 30 derniers jours)
Ali Almakhmari le 13 Jan 2023
Commenté : Ali Almakhmari le 13 Jan 2023
Hello everyone. I have a problem that I hope I can get some help on. So I have the following set of differential equations:
I need to solve for θ, ϕ, and ψ from those three differential equations. To be specific I need to solve for a value (not a function) of θ, ϕ, and ψ at a given moment. I have a for loop that runs, and at each iteration the value of p, q, and r change, so at each iteration I will need to do this integration and solve for a value of θ, ϕ, and ψ. What information do I have avaliable for me? The value of p, q, and r at each iteration, and I also know the initial values (iteration = 1) of θ, ϕ, and ψ, and finally, I know the time at each iteration as well. Is this possible in MATLAB? I tried implementing this with ODE45 and struggled very much. Any help is super appreciated.
##### 6 commentairesAfficher 4 commentaires plus anciensMasquer 4 commentaires plus anciens
Jan le 13 Jan 2023
Modifié(e) : Jan le 13 Jan 2023
It is still not clear, what you call "iteration". ODE45 solves the ODE in steps, so this is an iterative process. But steps can be rejected or re-calculated, such that "changes in each iteration" are not meaningful. If p,q,r are functions of time, there is no problem.
Please post the existing code, even if it fails. Then it might get clear immediately, how p,q,r are defined.
Ali Almakhmari le 13 Jan 2023
I tried so many things that I keep changing. This is the code I have now:
inc = 0.1;
time = 0:inc:100;
% the initial condition for the phi, theta, and psi
phi_sol = 0;
theta_sol = 30;
psi_sol = 0;
%for loop begins
for iter = 1:1:length(time)
omega_B(iter,:) = [0.01*time(iter), 0.1*time(iter), 0.001*time(iter)]; % [p, q, r]
syms thetas(t) phis(t) psis(t)
phi_dot = diff(phis) == omega_B(iter,1) + tan(thetas)*sin(phis)*omega_B(iter,2) + tan(thetas)*cos(phis)*omega_B(iter,3);
theta_dot = diff(thetas) == cos(phis)*omega_B(iter,2) - sin(phis)*omega_B(iter,3);
psi_dot = diff(psis) == sec(thetas)*sin(phis)*omega_B(iter,2) + sec(thetas)*cos(phis)*omega_B(iter,3);
odes = [phi_dot; theta_dot; psi_dot];
[phi_sol(iter+1), theta_sol(iter+1), psi_sol(iter+1)] = dsolve(odes, [phi_sol(iter), theta_sol(iter), psi_sol(iter)])
end
Note that omega_B represents p,q,r and I put random changing numbers for them instead of the real code because the real code is so big and complicated.

Connectez-vous pour commenter.

### Réponse acceptée

Steven Lord le 13 Jan 2023
The problem you've written out looks quite similar to the "Pass Extra Parameters to ODE Function" example on the ode45 documentation page. You have a system of three first-order ODEs rather than the system of two first-order ODEs given in the example and you have three extra parameters (p, q, and r) rather than the two (A and B) from the example, but those will require only a small modification to the odefun and to the ode45 call.
The fact that p, q, and r are functions of time rather than fixed values is a little bit more challenging, but the "ODE with Time-Dependent Terms" example on that same documentation page shows how you can handle them if you have a vector of data. If you have function handles for p, q, and r that's even easier.
f = @(t) t.^2;
[t, y] = ode45(@(t, y) myode(t, y, f), [0 1], 1);
plot(t, y)
Let's check the symbolic solution.
syms ts ys(ts) % Using ts and ys instead of t and y so I can use t and y for plotting below
solution = dsolve(diff(ys, ts) == ts^2*ys(ts), ys(0)==1)
solution =
And now plot the numeric solution and the value of the symbolic solution at the same times.
figure
plot(t, y, '+', t, subs(solution, ts, t), 'o')
Those look to be in pretty good agreement.
function dydt = myode(t, y, f)
fvalue = f(t);
dydt = fvalue*y; % dydt = t^2*y
end
##### 0 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

### Catégories

En savoir plus sur Programming dans Help Center et File Exchange

R2021a

### Community Treasure Hunt

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

Start Hunting!

Translated by