Heaviside function in ODE

Question

Pascal Schulthess 2018-4-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/394121-heaviside-function-in-ode

评论： Sayan Batabyal 2020-12-14

I have the following set of ODEs

dxdt(1) = -ka*x(1);
dxdt(2) = ka/V*x(1) - ke*x(2);
dxdt(3) = -a1*exp(-b1*x(3)).*x(3) + (x(2) - c1)*H(x(2) - c1);

wherein H(x) = 1 for x >= 0 and 0 otherwise is a Heaviside function. I've learned that such an ODE should be solved by using events.

This is what I currently have:

function mcm
% Time
tstart = 0;
tfinal = 2*24;
tspan = tstart:0.1:tfinal;
% Initial conditions
x0 = [200; 0; 0];
% Parameters
ka = 2.4;
V = 14;
ke = 0.39;
a1 = 0.7;
b1 = 0.31;
c1 = 3.7;
% Options with event function
options = odeset('Events', @heavi);
% Solve until first terminal event
[t, x, te, xe, ie] = ode45(@model, tspan, x0, options);
% Plot
plotnames = ['x1'; 'x2'; 'x3'];
for i = 1:length(x0)
    subplot(3, 1, i);
    plot(t, x(:, i))
    title(plotnames(i, :));
end
% Model
    function dxdt = model(t, x)
          % Initialise
          dxdt = zeros(3, 1);
          % Model equations
          dxdt(1) = -ka*x(1);
          dxdt(2) = ka/V*x(1) - ke*x(2);
          dxdt(3) = -a1*exp(-b1*x(3)).*x(3);
      end
% Event function
    function [value, isterminal, direction] = heavi(t, x)
        value = x(2) - c1;
        isterminal = 0;
        direction = 0;
    end
end

But I'm not sure how to add the Heaviside part to the third ODE. Any help would be greatly appreciated.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2018-4-11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/394121-heaviside-function-in-ode#answer_314539

1. Set isterminal=1 in the event function and organize your code as in the "ballode" example.

2. Switch a variable "flag" in "mcm" from 0 to 1 and vice versa each time you return from ode45 (i.e. when the event function detects a zero crossing of "value") and set

dxdt(3) = -a1*exp(-b1*x(3)).*x(3)+flag*(x(2)-c1);

Best wishes

Torsten.