MATLAB Answers

Heaviside function in ODE

23 views (last 30 days)
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.

Accepted Answer

Torsten
Torsten on 11 Apr 2018
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.
  2 Comments
Sayan Batabyal
Sayan Batabyal on 14 Dec 2020
@Torsten, Do you know how to change the value of flag inside the events function? If you can, please elaboarte with a small code snippet.

Sign in to comment.

More Answers (0)

Tags

Products

Community Treasure Hunt

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

Start Hunting!

Translated by