Time-based integration of ODEs

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.

[Davide Masiello]

This seems like a very standard problem and could be solved easily in MatLab.

I suggest you post your code so that we can suggest what to change/improve.

[Star Strider]

What actually defines an ‘iteration’?

Do the ‘tspan’ limits change in each iteration? If so, how are they defined?

Are p, q, and r defined as function results (and if so, what is the argument to the function) or are they elements of vectors?

[Jan]

@Ali Almakhmari: What does "each iteration" mean? Are p,q,r functions of the time? If not, how are they defined?

[Ali Almakhmari]

Wow. Thank you all for your questions. And yes p, q, and r are a function of time and they change in each iteration. Each iteration represents an increment in time. While tspan can be defined in any way by the user before entering the for loop.

[Jan]

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]

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.

Answer 1

[Steven Lord]

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