Heaviside function in ODE
7 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Pascal Schulthess
am 11 Apr. 2018
Kommentiert: Sayan Batabyal
am 14 Dez. 2020
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 Kommentare
Akzeptierte Antwort
Torsten
am 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 Kommentare
Arjun Chakrawal
am 6 Feb. 2020
Pascal Schulthess-- could you please share your code on how did you actually implemented the suggestions by Torsten?
Sayan Batabyal
am 14 Dez. 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.
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Produkte
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!