Skip to content

Commit be1079b

Browse files
Steven-RobertsSteven-Roberts
authored
Feature/oregonator docs (#119)
* Start Oregonator docs * Add Canonical constructor docs * Fix Jacobian
1 parent e4e4ebb commit be1079b

File tree

8 files changed

+242
-16
lines changed

8 files changed

+242
-16
lines changed

‎docs/references.bib‎

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -120,3 +120,38 @@ @article{PS19
120120
publisher={Copernicus Publications G{\"o}ttingen, Germany},
121121
doi={10.5194/npg-26-109-2019}
122122
}
123+
124+
@article{FN74,
125+
author = {Field, Richard J. and Noyes, Richard M.},
126+
title = "{Oscillations in chemical systems. {IV}. {L}imit cycle behavior in a model of a real chemical reaction}",
127+
journal = {The Journal of Chemical Physics},
128+
volume = {60},
129+
number = {5},
130+
pages = {1877-1884},
131+
year = {1974},
132+
issn = {0021-9606},
133+
doi = {10.1063/1.1681288}
134+
}
135+
136+
@incollection{EH76,
137+
title = {Comparing Numerical Methods for the Solution of Stiff Systems of {ODE}s Arising in Chemistry},
138+
editor = {L. LAPIDUS and W.E. SCHIESSER},
139+
booktitle = {Numerical Methods for Differential Systems},
140+
publisher = {Academic Press},
141+
pages = {45-66},
142+
year = {1976},
143+
isbn = {978-0-12-436640-4},
144+
doi = {10.1016/B978-0-12-436640-4.50008-3},
145+
author = {W.H. Enright and T.E. Hull}
146+
}
147+
148+
@article{GW82,
149+
author = {Björn A. Gottwald and Gerhard Wanner},
150+
title ={Comparison of Numerical Methods for Stiff Differential Equations in Biology and Chemistry},
151+
journal = {SIMULATION},
152+
volume = {38},
153+
number = {2},
154+
pages = {61-66},
155+
year = {1982},
156+
doi = {10.1177/003754978203800206}
157+
}

‎toolbox/+otp/+oregonator/+presets/Canonical.m‎

Lines changed: 40 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,48 @@
11
classdef Canonical < otp.oregonator.OregonatorProblem
2+
% The Oregonator configuration used in :cite:p:`HW96` (p. 144) called OREGO. It uses time span $t \in [0, 360]$,
3+
% initial condition $y_0 = [1, 2, 3]^T$, and parameters
4+
%
5+
% $$
6+
% f &= 1, \\
7+
% q &= 8.375 \cdot 10^{-6}, \\
8+
% s &= 77.27, \\
9+
% w &= 0.1610.
10+
% $$
11+
212
methods
3-
function obj = Canonical
13+
function obj = Canonical(varargin)
14+
% Create the Canonical problem object.
15+
%
16+
% Parameters
17+
% ----------
18+
% varargin
19+
% A variable number of name-value pairs. The accepted names are
20+
%
21+
% - ``f`` – Value of $f$.
22+
% - ``q`` – Value of $q$.
23+
% - ``s`` – Value of $s$.
24+
% - ``w`` – Value of $w$.
25+
%
26+
% Returns
27+
% -------
28+
% obj : Canonical
29+
% The constructed problem.
30+
31+
p = inputParser;
32+
p.addParameter('f', 1);
33+
p.addParameter('q', 8.375e-6);
34+
p.addParameter('s', 77.27);
35+
p.addParameter('w', 0.1610);
36+
p.parse(varargin{:});
37+
opts = p.Results;
38+
439
tspan = [0, 360];
540
y0 = [1; 2; 3];
641
params = otp.oregonator.OregonatorParameters;
42+
params.F = 1;
43+
params.Q = 8.375e-6;
44+
params.S = 77.27;
45+
params.W = 0.1610;
746

847
obj = obj@otp.oregonator.OregonatorProblem(tspan, y0, params);
948
end

‎toolbox/+otp/+oregonator/+presets/EnrightHull.m‎

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,35 @@
1+
classdef EnrightHull < otp.oregonator.OregonatorProblem
2+
% The Oregonator configuration used in problem CHM 9 of :cite:p:`EH76`. It uses time span $t \in [0, 300]$, initial
3+
% condition $y_0 = [4, 1.1, 4]^T$, and parameters
4+
%
5+
% $$
6+
% f &= 1, \\
7+
% q &= 8.375 \cdot 10^{-6}, \\
8+
% s &= 77.27, \\
9+
% w &= 0.1610.
10+
% $$
11+
12+
methods
13+
function obj = EnrightHull
14+
% Create the EnrightHull problem object.
15+
%
16+
% Parameters
17+
% ----------
18+
%
19+
% Returns
20+
% -------
21+
% obj : EnrightHull
22+
% The constructed problem.
23+
24+
tspan = [0, 300];
25+
y0 = [4; 1.1; 4];
26+
params = otp.oregonator.OregonatorParameters;
27+
params.F = 1;
28+
params.Q = 8.375e-6;
29+
params.S = 77.27;
30+
params.W = 0.1610;
31+
32+
obj = obj@otp.oregonator.OregonatorProblem(tspan, y0, params);
33+
end
34+
end
35+
end

‎toolbox/+otp/+oregonator/+presets/GottwaldWanner.m‎

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
classdef GottwaldWanner < otp.oregonator.OregonatorProblem
2+
% The Oregonator configuration from :cite:p:`GW82` which uses time span $t \in [0, 302.85805]$, initial condition
3+
% $y_0 = [4, 1.331391, 2.852348]^T$, and parameters
4+
%
5+
% $$
6+
% f &= 1, \\
7+
% q &= 8.375 \cdot 10^{-6}, \\
8+
% s &= 77.27, \\
9+
% w &= 0.1610.
10+
% $$
11+
%
12+
% The initial condition lies on a stable limit cycle, and the time span is chosen to be one period.
13+
14+
methods
15+
function obj = GottwaldWanner
16+
% Create the GottwaldWanner problem object.
17+
%
18+
% Parameters
19+
% ----------
20+
%
21+
% Returns
22+
% -------
23+
% obj : GottwaldWanner
24+
% The constructed problem.
25+
26+
tspan = [0, 302.85805];
27+
y0 = [4; 1.331391; 2.852348];
28+
params = otp.oregonator.OregonatorParameters;
29+
params.F = 1;
30+
params.Q = 8.375e-6;
31+
params.S = 77.27;
32+
params.W = 0.1610;
33+
34+
obj = obj@otp.oregonator.OregonatorProblem(tspan, y0, params);
35+
end
36+
end
37+
end

‎toolbox/+otp/+oregonator/OregonatorParameters.m‎

Lines changed: 14 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,17 @@
11
classdef OregonatorParameters
2-
%OREGONATORPARAMETERS
2+
% Parameters for the Oregonator problem.
3+
properties
4+
% The stoichiometric factor $f$.
5+
F %MATLAB ONLY: (1, 1) {mustBeReal, mustBeFinite}
6+
7+
% The reaction constant $q$.
8+
Q %MATLAB ONLY: (1, 1) {mustBeReal, mustBeFinite}
9+
10+
% The reaction constant $s$.
11+
S %MATLAB ONLY: (1, 1) {mustBeReal, mustBeFinite}
12+
13+
% The reaction constant $w$.
14+
W %MATLAB ONLY: (1, 1) {mustBeReal, mustBeFinite}
15+
end
316
end
417

‎toolbox/+otp/+oregonator/OregonatorProblem.m‎

Lines changed: 73 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,25 +1,92 @@
11
classdef OregonatorProblem < otp.Problem
2+
% A periodic, three-variable model for the Belousov–Zhabotinsky reaction.
3+
%
4+
% The Oregonator chemical reaction :cite:p:`FN74` is given by
5+
%
6+
%
7+
% $$
8+
% \ce{
9+
% A + Y &-> X + P \\
10+
% X + Y &-> 2P \\
11+
% A + X &-> 2X + 2Z \\
12+
% 2X &-> A + P \\
13+
% B + Z &-> 1/2 $f$ Y.
14+
% }
15+
% $$
16+
%
17+
% Species $\ce{A = BrO3-}$, $\ce{B = CH2(COOH)2}$, and $\ce{P = HOBr}$ or $\ce{BrCH(COOH)2}$ are assumed to be
18+
% constant. The dynamic concentrations of intermediate species $\ce{X = HBrO2}$, $\ce{Y = Br-}$, and
19+
% $\ce{Z = Ce(IV)}$ can be modeled by an ODE in three variables. Field and Noyes :cite:p:`FN74` proposed the
20+
% nondimensionalized form
21+
%
22+
% $$
23+
% α' &= s(η - η α + α - q α^2), \\
24+
% η' &= s^{-1}(-η - η α + f ρ), \\
25+
% ρ' &= w (α - ρ),
26+
% $$
27+
%
28+
% where $α$, $η$, and $ρ$ are scaled concentrations of $\ce{X}$, $\ce{Y}$, and $\ce{Z}$, respectively.
29+
%
30+
% Notes
31+
% -----
32+
% +---------------------+------------------------------------------------+
33+
% | Type | ODE |
34+
% +---------------------+------------------------------------------------+
35+
% | Number of Variables | 3 |
36+
% +---------------------+------------------------------------------------+
37+
% | Stiff | typically, depending on $f$, $q$, $s$, and $w$ |
38+
% +---------------------+------------------------------------------------+
39+
%
40+
% Example
41+
% -------
42+
% >>> problem = otp.oregonator.presets.Canonical;
43+
% >>> sol = problem.solve();
44+
% >>> problem.plotPhaseSpace(sol, 'Vars', [2, 3]);
45+
246
methods
347
function obj = OregonatorProblem(timeSpan, y0, parameters)
48+
% Create a Oregonator problem object.
49+
%
50+
% Parameters
51+
% ----------
52+
% timeSpan : numeric(1, 2)
53+
% The start and final time.
54+
% y0 : numeric(3, 1)
55+
% The initial conditions.
56+
% parameters : OregonatorParameters
57+
% The parameters.
58+
%
59+
% Returns
60+
% -------
61+
% obj : OregonatorProblem
62+
% The constructed problem.
463
obj@otp.Problem('Oregonator', 3, timeSpan, y0, parameters);
564
end
665
end
766

867
methods (Access = protected)
968
function onSettingsChanged(obj)
10-
obj.RHS = otp.RHS(@otp.oregonator.f, ...
11-
'Jacobian', @otp.oregonator.jacobian, ...
69+
f = obj.Parameters.F;
70+
q = obj.Parameters.Q;
71+
s = obj.Parameters.S;
72+
w = obj.Parameters.W;
73+
74+
obj.RHS = otp.RHS(@(t, y) otp.oregonator.f(t, y, f, q, s, w), ...
75+
'Jacobian', @(t, y) otp.oregonator.jacobian(t, y, f, q, s, w), ...
1276
'Vectorized', 'on');
1377
end
1478

1579
function fig = internalPlot(obj, t, y, varargin)
16-
fig = internalPlot@otp.Problem(obj, t, y, ...
17-
'yscale', 'log', varargin{:});
80+
fig = internalPlot@otp.Problem(obj, t, y, 'yscale', 'log', varargin{:});
81+
end
82+
83+
function fig = internalPlotPhaseSpace(obj, t, y, varargin)
84+
fig = internalPlotPhaseSpace@otp.Problem(obj, t, y, 'xscale', 'log', 'yscale', 'log', 'zscale', 'log', ...
85+
varargin{:});
1886
end
1987

2088
function mov = internalMovie(obj, t, y, varargin)
21-
mov = internalMovie@otp.Problem(obj, t, y, ...
22-
'yscale', 'log', varargin{:});
89+
mov = internalMovie@otp.Problem(obj, t, y, 'yscale', 'log', varargin{:});
2390
end
2491
end
2592
end

‎toolbox/+otp/+oregonator/f.m‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,12 @@
1-
function dy = f(~, y)
1+
function dy = f(~, y, f, q, s, w)
22

33
y1 = y(1, :);
44
y2 = y(2, :);
55
y3 = y(3, :);
66

77
dy = [ ...
8-
77.27 * (y2 + y1 .* (1 - 8.375e-6 * y1 - y2)); ...
9-
(y3 - (1 + y1) .* y2) / 77.27; ...
10-
0.161 * (y1 - y3)];
8+
s * (y2 + y1 .* (1 - q * y1 - y2)); ...
9+
(f * y3 - (1 + y1) .* y2) / s; ...
10+
w * (y1 - y3)];
1111

1212
end

‎toolbox/+otp/+oregonator/jacobian.m‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,11 @@
1-
function J = jacobian(~, y)
1+
function J = jacobian(~, y, f, q, s, w)
22

33
y1 = y(1);
44
y2 = y(2);
55

66
J = [ ...
7-
77.27 * (1 - y2 - 1.675e-5 * y1), 77.27 * (1 - y1), 0; ...
8-
-y2 / 77.27, -(1 + y1) / 77.27, 1 / 77.27; ...
9-
0.161, 0, -0.161];
7+
s * (1 - 2 * q * y1 - y2), s * (1 - y1), 0; ...
8+
-y2 / s, -(1 + y1) / s, f / s; ...
9+
w, 0, -w];
1010

1111
end

0 commit comments

Comments
 (0)