MATLAB Answers

Heaviside function in ODE

31 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.

  0 Comments

Sign in to comment.

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.

  1 Comment

Arjun Chakrawal
Arjun Chakrawal on 6 Feb 2020
Pascal Schulthess-- could you please share your code on how did you actually implemented the suggestions by Torsten?

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!