How to call an ode solver within another function?

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?

[Walter Roberson]

[~,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.]

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?

Answer 1

[Sam Chak]

Hi @L'O.G.

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