How to call an ode solver within another function?

조회 수: 6 (최근 30일)
L'O.G.
L'O.G. 2023년 10월 5일
답변: Sam Chak 2023년 10월 5일
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
Walter Roberson
Walter Roberson 2023년 10월 5일
[~,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. 2023년 10월 5일
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?

댓글을 달려면 로그인하십시오.

채택된 답변

Sam Chak
Sam Chak 2023년 10월 5일
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개)

카테고리

Help CenterFile Exchange에서 Programming에 대해 자세히 알아보기

제품


릴리스

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by