MATLAB Answers

0

Heaviside function in ODE

Asked by Pascal Schulthess on 11 Apr 2018
Latest activity Answered by Torsten
on 11 Apr 2018
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.

Tags

Products

1 Answer

Torsten
Answer by Torsten
on 11 Apr 2018
 Accepted Answer

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.

  0 Comments

Sign in to comment.