Time-based integration of ODEs

10 次查看（过去 30 天）
Ali Almakhmari
Ali Almakhmari 2023-1-13
评论： Ali Almakhmari 2023-1-13
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 个评论
显示 4更早的评论
Jan
Jan 2023-1-13
编辑：Jan 2023-1-13
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
Ali Almakhmari 2023-1-13
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Steven Lord
Steven Lord 2023-1-13
Ran in:
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 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by