How to call an ode solver within another function?

6 Ansichten (letzte 30 Tage)
L'O.G. am 5 Okt. 2023
Beantwortet: Sam Chak am 5 Okt. 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
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 Kommentare
Walter Roberson
Walter Roberson am 5 Okt. 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. am 5 Okt. 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?

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Sam Chak
Sam Chak am 5 Okt. 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

Weitere Antworten (0)


Mehr zu Programming finden Sie in Help Center und File Exchange




Community Treasure Hunt

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

Start Hunting!

Translated by