Clear Filters
Clear Filters

How to call an ode solver within another function?

9 views (last 30 days)
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?
Walter Roberson
Walter Roberson on 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. on 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?

Sign in to comment.

Accepted Answer

Sam Chak
Sam Chak on 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

More Answers (0)


Find more on Programming in Help Center and File Exchange




Community Treasure Hunt

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

Start Hunting!