How to call an ode solver within another function?

6 vues (au cours des 30 derniers jours)
L'O.G.
L'O.G. le 5 Oct 2023
I solve dx/dt = x(t) using an ode solver
sol = ode45(@(t,x) u(t,x,ode,eps,a0),[t0 tf],[x1;x2]);
that calls a function u for the velocity
function v = u(t,x,ode,eps,a0)
deltaX = x2 - x1; % they are different at each time point: x2(t) and x1(t)
[~,at] = ode45(@(t,a_t) ode(t,a_t,deltaX),[0 t + eps],a0);
a_t = interp1(linspace(0, t + eps, numel(at)), at, t, 'linear');
% ... some calculation involving a_t
end
where I wish to solve the differential equation:
ode = @(t,a_t,deltaX) a0 * (x2-x1) / a_t);
which gives da/dt for a(t) for each time point since x2 and x1 are updated at each time point.
Since I get an error message saying the last entry in tspan must be different from the first entry, I introduced a small positive number eps and then interpolate for the given t, but this interpolated value a_t changes depending on the tspan value, so I don't believe I did this correctly.
In other words, I am trying to solve for x(t) given the velocity that itself depends on a term a(t) that depends on deltaX = x2(t)-x1(t). I feel like the solver for a(t) needs to be within the function u since a(t) is required for other calculations at a given time. How do I do this?
  2 commentaires
Walter Roberson
Walter Roberson le 5 Oct 2023
[~,at] = ode45(@(t,a_t) ode(t,a_t,deltaX),[0 t + eps],a0);
a_t = interp1(linspace(0, t + eps, numel(at)), at, t, 'linear');
That code is wrong. Use
tspan = linspace(0, t+eps, 5);
sol = ode45(@(t,a_t) ode(t,a_t,deltaX), tspan, a0);
a_t = deval(sol, t);
L'O.G.
L'O.G. le 5 Oct 2023
Thanks. I get the same result. Does it make sense to have the ode solver for a(t) nested in the other function, and for tspan to go from 0 to t+eps?

Connectez-vous pour commenter.

Réponse acceptée

Sam Chak
Sam Chak le 5 Oct 2023
I'm uncertain if I've correctly deduced the mathematical problem from your code. However, a singularity will occur if due to a division-by-zero term in . Since these are coupled ordinary differential equations (ODEs), if the integrations of occur at the same time interval, you can include them in the same odefcn() function, as demonstrated below:
t0 = 0;
tf = 400;
tspan = [t0 tf];
a0 = 0.1;
x0 = [0 1 a0];
[t, x] = ode45(@(t, x) odefcn(t, x, a0), tspan, x0);
plot(t, x), grid on
xlabel({'$t$'}, 'interpreter', 'latex', fontsize=14),
ylabel({'$\mathbf{x}(t)$'}, 'interpreter', 'latex', fontsize=14)
legend({'$x_{1}(t)$', '$x_{2}(t)$', '$a_{t}$'}, 'interpreter', 'latex', fontsize=16)
function dxdt = odefcn(t, x, a0)
dxdt = zeros(3, 1);
at = x(3);
xref = 0.25*sin(pi/5*t) + 0.1;
dxdt(1) = x(2) - x(1); % sample ode 1
dxdt(2) = - 0.1*(x(1) - xref) - 0.2*at*tanh((x(2) - x(1))/0.01); % sample ode 2
dxdt(3) = a0*(x(2) - x(1))/at; % singularity happens when at = 0 % user-supplied ode
end

Plus de réponses (0)

Produits


Version

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by