How to call an ode solver within another function?
6 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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
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);
Réponse acceptée
Sam Chak
le 5 Oct 2023
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
0 commentaires
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!