Skip to content

Commit 2bcaf67

Browse files
Steven-RobertsSteven-Roberts
authored
Docs/hires (#135)
* Add parameters to HIRES * Finish HIRES docs
1 parent 5b943b3 commit 2bcaf67

File tree

8 files changed

+274
-81
lines changed

8 files changed

+274
-81
lines changed

‎docs/references.bib

+33-3
Original file line numberDiff line numberDiff line change
@@ -173,9 +173,39 @@ @article{Pet86
173173
number = {4},
174174
pages = {837--852},
175175
publisher = {Society for Industrial and Applied Mathematics},
176-
title = {Order Results for Implicit Runge-Kutta Methods Applied to Differential/ Algebraic Systems},
176+
title = {Order Results for Implicit {Runge--Kutta} Methods Applied to Differential/Algebraic Systems},
177177
volume = {23},
178178
year = {1986},
179-
url = {https://www.jstor.org/stable/2157625}
179+
doi = {10.1137/0723054}
180+
}
181+
182+
@article{Sch75,
183+
title = {A new approach to explain the ``high irradiance responses'' of photomorphogenesis on the basis of phytochrome},
184+
volume = {2},
185+
ISSN = {1432-1416},
186+
DOI = {10.1007/bf00276015},
187+
number = {1},
188+
journal = {Journal of Mathematical Biology},
189+
author = {Schäfer, E.},
190+
year = {1975},
191+
pages = {41-56}
192+
}
193+
194+
@article{Got77,
195+
title = {{MISS} - ein einfaches Simulations-System für biologische und chemische Prozesse},
196+
volume = {3},
197+
journal = {EDV in Medizin und Biologie},
198+
author = {B. A. Gottwald},
199+
year = {1977},
200+
pages = {85-90}
201+
}
202+
203+
@article{dSL98,
204+
title={Collecting real-life problems to test solvers for implicit differential equations},
205+
author={de Swart, Jacques J. B. and Lioen, Walter M.},
206+
journal={CWI Quarterly},
207+
volume={11},
208+
number={1},
209+
pages={83-100},
210+
year={1998}
180211
}
181-

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

+33-1
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,41 @@
11
classdef Canonical < otp.hires.HIRESProblem
2+
% HIRES preset proposed in :cite:p:`HW96` (pp. 144-145) which uses time span $t ∈ [0, 321.8122]$, initial conditions
3+
%
4+
% $$
5+
% y(0) = [1, 0, 0, 0, 0, 0, 0, 0.0057]^T,
6+
% $$
7+
%
8+
% and parameters
9+
%
10+
% $$
11+
% k_1 &= 1.71, \quad & k_2 &= 0.43, \quad & k_3 &= 8.32, \quad & k_4 &= 0.69, \quad & k_5 &= 0.035, \\
12+
% k_6 &= 8.32, & k_{+} &= 280, & k_{-} &= 0.69, & k^* &= 0.69, & o_{k_s} &= 0.0007.
13+
% $$
14+
215
methods
316
function obj = Canonical(varargin)
17+
% Create the canonical HIRES problem object.
18+
%
19+
% Parameters
20+
% ----------
21+
% varargin
22+
% A variable number of name-value pairs. The accepted names are
23+
%
24+
% - ``K1`` – Value of $k_1$.
25+
% - ``K2`` – Value of $k_2$.
26+
% - ``K3`` – Value of $k_3$.
27+
% - ``K4`` – Value of $k_4$.
28+
% - ``K5`` – Value of $k_5$.
29+
% - ``K6`` – Value of $k_6$.
30+
% - ``KPlus`` – Value of $k_{+}$.
31+
% - ``KMinus`` – Value of $k_{-}$.
32+
% - ``KStar`` – Value of $k^*$.
33+
% - ``OKS`` – Value of $o_{k_s}$.
34+
435
tspan = [0, 321.8122];
536
y0 = [1; 0; 0; 0; 0; 0; 0; 0.0057];
6-
params = otp.hires.HIRESParameters(varargin{:});
37+
params = otp.hires.HIRESParameters('K1', 1.71, 'K2', 0.43, 'K3', 8.32, 'K4', 0.69, 'K5', 0.035, ...
38+
'K6', 8.32, 'KPlus', 280, 'KMinus', 0.69, 'KStar', 0.69, 'OKS', 0.0007, varargin{:});
739
obj = obj@otp.hires.HIRESProblem(tspan, y0, params);
840
end
941
end

‎toolbox/+otp/+hires/HIRESParameters.m

+32
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,38 @@
11
classdef HIRESParameters < otp.Parameters
22
% Parameters for the HIRES problem.
33

4+
properties
5+
% The reaction rate $k_1$.
6+
K1 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical}
7+
8+
% The reaction rate $k_2$.
9+
K2 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical}
10+
11+
% The reaction rate $k_3$.
12+
K3 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical}
13+
14+
% The reaction rate $k_4$.
15+
K4 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical}
16+
17+
% The reaction rate $k_5$.
18+
K5 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical}
19+
20+
% The reaction rate $k_6$.
21+
K6 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical}
22+
23+
% The reaction rate $k_{+}$.
24+
KPlus %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical}
25+
26+
% The reaction rate $k_{-}$.
27+
KMinus %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical}
28+
29+
% The reaction rate $k^*$.
30+
KStar %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical}
31+
32+
% The source term $o_{k_s}$.
33+
OKS %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical}
34+
end
35+
436
methods
537
function obj = HIRESParameters(varargin)
638
% Create a HIRES parameters object.

‎toolbox/+otp/+hires/HIRESProblem.m

+104-5
Original file line numberDiff line numberDiff line change
@@ -1,18 +1,117 @@
1-
21
classdef HIRESProblem < otp.Problem
2+
% An eight-variable high irradiance response model from plant physiology.
3+
%
4+
% In :cite:p:`Sch75`, the following chemical reactions were proposed to explain the high irradiance responses of
5+
% photomorphogenesis on the basis of phytochrome:
6+
%
7+
% $$
8+
% \begin{array}{cccc}
9+
% \xrightarrow{o_{k_s}} & \ce{P_{{r}}} & \xrightleftharpoons[k_2]{k_1} & \ce{P_{fr}} \\
10+
% & {\scriptsize k_6} ↑ & & ↓ {\scriptsize k_3} \\
11+
% & \ce{P_{{r}}X} & \xrightleftharpoons[k_2]{k_1} & \ce{P_{fr}X} \\
12+
% & {\scriptsize k_5} ↑ & & ↓ {\scriptsize k_4} \\
13+
% & \ce{P_{{r}}X'} & \xrightleftharpoons[k_2]{k_1} & \ce{P_{fr}X'} \\
14+
% \end{array}
15+
% $$
16+
%
17+
% $$
18+
% \begin{array}{ccccc}
19+
% \ce{E + P_{{r}}X'} & \xleftarrow{k_2} & \ce{P_{fr}X'E} & \xrightleftharpoons[k_{-}]{k_{+}} & \ce{P_{fr}X' + E} \\
20+
% & & ↓ {\scriptsize k^*} \\
21+
% & & \ce{P_{fr}' + X' + E}
22+
% \end{array}
23+
% $$
24+
%
25+
% Reactants $\ce{P_{{r}}}$ and $\ce{P_{fr}}$ are the red and far-red absorbing form of phytochrome, respectively.
26+
% These can be bound by receptors $\ce{X}$ and $\ce{X'}$, partially influenced by enzyme $\ce{E}$. The system is
27+
% modeled by the differential equations :cite:p:`Got77` :cite:p:`dSL98`
28+
%
29+
% $$
30+
% \frac{d}{dt} \ce{[P_{{r}}]} &= -k_1 \ce{[P_{{r}}]} + k_2 \ce{[P_{fr}]} + k_6 \ce{[P_{{r}}X]} + o_{k_s}, \\
31+
% \frac{d}{dt} \ce{[P_{fr}]} &= k_1 \ce{[P_{{r}}]} - (k_2 + k_3) \ce{[P_{fr}]}, \\
32+
% \frac{d}{dt} \ce{[P_{{r}}X]} &= -(k_1 + k_6) \ce{[P_{{r}}X]} + k_2 \ce{[P_{fr}X]} + k_5 \ce{[P_{fr}X']}, \\
33+
% \frac{d}{dt} \ce{[P_{fr}X]} &= k_3 \ce{[P_{fr}]} + k_1 \ce{[P_{{r}}X]} - (k_2 + k_4) \ce{[P_{fr}X]}, \\
34+
% \frac{d}{dt} \ce{[P_{{r}}X']} &= -(k_1 + k_5) \ce{[P_{{r}}X']} + k_2 \ce{[P_{fr}X']} + k_2 \ce{[P_{fr}X'E]}, \\
35+
% \frac{d}{dt} \ce{[P_{fr}X']} &= k_4 \ce{[P_{fr}X]} + k_1 \ce{[P_{{r}}X']} - k_2 \ce{[P_{fr}X']}
36+
% + k_{-} \ce{[P_{fr}X'E]} - k_{+} \ce{[P_{fr}X']} \ce{[E]}, \\
37+
% \frac{d}{dt} \ce{[P_{fr}X'E]} &= -(k_2 + k_{-} + k^*) \ce{[P_{fr}X'E]} + k_{+} \ce{[P_{fr}X']} \ce{[E]}, \\
38+
% \frac{d}{dt} \ce{[E]} &= (k_2 + k_{-} + k^*) \ce{[P_{fr}X'E]} - k_{+} \ce{[P_{fr}X']} \ce{[E]}.
39+
% $$
40+
%
41+
% Notes
42+
% -----
43+
% +---------------------+--------------------------------------------------------------------+
44+
% | Type | ODE |
45+
% +---------------------+--------------------------------------------------------------------+
46+
% | Number of Variables | 8 |
47+
% +---------------------+--------------------------------------------------------------------+
48+
% | Stiff | typically, depending on $k_1, …, k_6$, $k_{+}$, $k_{-}$, and $k^*$ |
49+
% +---------------------+--------------------------------------------------------------------+
50+
%
51+
% Example
52+
% -------
53+
% >>> problem = otp.hires.presets.Canonical;
54+
% >>> sol = problem.solve('AbsTol', 1e-12);
55+
% >>> problem.plot(sol);
356

457
methods
558
function obj = HIRESProblem(timeSpan, y0, parameters)
59+
% Create a HIRES problem object.
60+
%
61+
% Parameters
62+
% ----------
63+
% timeSpan : numeric(1, 2)
64+
% The start and final time.
65+
% y0 : numeric(2, 1)
66+
% The initial conditions.
67+
% parameters : otp.hires.HIRESParameters
68+
% The parameters.
669
obj@otp.Problem('HIRES', 8, timeSpan, y0, parameters);
770
end
871
end
972

1073
methods (Access = protected)
1174
function onSettingsChanged(obj)
12-
obj.RHS = otp.RHS(@otp.hires.f, ...
13-
'Jacobian', @otp.hires.jacobian, ...
14-
'JacobianVectorProduct', @otp.hires.jacobianVectorProduct, ...
15-
'JacobianAdjointVectorProduct', @otp.hires.jacobianAdjointVectorProduct);
75+
k1 = obj.Parameters.K1;
76+
k2 = obj.Parameters.K2;
77+
k3 = obj.Parameters.K3;
78+
k4 = obj.Parameters.K4;
79+
k5 = obj.Parameters.K5;
80+
k6 = obj.Parameters.K6;
81+
kPlus = obj.Parameters.KPlus;
82+
kMinus = obj.Parameters.KMinus;
83+
kStar = obj.Parameters.KStar;
84+
oks = obj.Parameters.OKS;
85+
86+
obj.RHS = otp.RHS(@(t, y) otp.hires.f(t, y, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, oks), ...
87+
'Jacobian', @(t, y) otp.hires.jacobian(t, y, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, oks), ...
88+
'JacobianVectorProduct', @(t, y, v) ...
89+
otp.hires.jacobianVectorProduct(t, y, v, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, oks), ...
90+
'JacobianAdjointVectorProduct', @(t, y, v) ...
91+
otp.hires.jacobianAdjointVectorProduct(t, y, v, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, oks), ...
92+
'Vectorized', 'on', ...
93+
'NonNegative', 1:obj.NumVars);
94+
end
95+
96+
function label = internalIndex2label(obj, index)
97+
switch index
98+
case 1
99+
label = 'P_r';
100+
case 2
101+
label = 'P_{fr}';
102+
case 3
103+
label = 'P_r X';
104+
case 4
105+
label = 'P_{fr} X';
106+
case 5
107+
label = 'P_r X''';
108+
case 6
109+
label = 'P_{fr} X''';
110+
case 7
111+
label = 'P_{fr} X'' E';
112+
case 8
113+
label = 'E';
114+
end
16115
end
17116

18117
function fig = internalPlot(obj, t, y, varargin)

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

+26-11
Original file line numberDiff line numberDiff line change
@@ -1,13 +1,28 @@
1-
function yPrime = f(~, y)
1+
function yPrime = f(~, y, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, oks)
22

3-
yPrime = [
4-
-1.71 * y(1) + 0.43 * y(2) + 8.32 * y(3) + 0.0007;
5-
1.71 * y(1) - 8.75 * y(2);
6-
-10.03 * y(3) + 0.43 * y(4) + 0.035 * y(5);
7-
8.32 * y(2) + 1.71 * y(3) - 1.12 * y(4);
8-
-1.745 * y(5) + 0.43 * y(6) + 0.43 * y(7);
9-
-280 * y(6) * y(8) + 0.69 * y(4) + 1.71 * y(5) - 0.43 * y(6) + 0.69 * y(7);
10-
280 * y(6) * y(8) - 1.81 * y(7);
11-
-280 * y(6) * y(8) + 1.81 * y(7)];
3+
k1Pr = k1 * y(1, :);
4+
k2Pfr = k2 * y(2, :);
5+
k3Pfr = k3 * y(2, :);
6+
k1PrX = k1 * y(3, :);
7+
k6PrX = k6 * y(3, :);
8+
k2PfrX = k2 * y(4, :);
9+
k4PfrX = k4 * y(4, :);
10+
k1PrXPrime = k1 * y(5, :);
11+
k5PrXPrime = k5 * y(5, :);
12+
k2PfrXPrime = k2 * y(6, :);
13+
k2PfrXPrimeE = k2 * y(7, :);
14+
kMinusPfrXPrimeE = kMinus * y(7, :);
15+
kPlusPfrXPrimeE = kPlus * y(6, :) * y(8, :);
16+
yPrime7 = -k2PfrXPrimeE - (kMinus + kStar) * y(7, :) + kPlusPfrXPrimeE;
1217

13-
end
18+
yPrime = [ ...
19+
-k1Pr + k2Pfr + k6PrX + oks;
20+
k1Pr - k2Pfr - k3Pfr;
21+
-k1PrX - k6PrX + k2PfrX + k5PrXPrime;
22+
k3Pfr + k1PrX - k2PfrX - k4PfrX;
23+
-k1PrXPrime - k5PrXPrime + k2PfrXPrime + k2PfrXPrimeE;
24+
k4PfrX + k1PrXPrime - k2PfrXPrime + kMinusPfrXPrimeE - kPlusPfrXPrimeE;
25+
yPrime7;
26+
-yPrime7];
27+
28+
end

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

+9-19
Original file line numberDiff line numberDiff line change
@@ -1,23 +1,13 @@
1-
function J = jacobian(~, y)
1+
function jac = jacobian(~, y, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, ~)
22

3-
y6 = y(6);
4-
y8 = y(8);
3+
kPlusPfrXPrime = kPlus * y(6);
4+
kPlusE = kPlus * y(8);
5+
kSum = k2 + kMinus + kStar;
56

6-
J = sparse( ...
7-
[1,2,1,2,4,1,3,4,3,4,6,3,5,6,5,6,7,8,5,6,7,8,6,7,8], ...
8-
[1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6,6,7,7,7,7,8,8,8], ...
9-
[-1.71, 1.71, 0.43, -8.75, 8.32, 8.32, -10.03, 1.71, 0.43, -1.12, 0.69, 0.035, -1.745, 1.71, ...
10-
0.43, -280 * y8 - 0.43, 280 * y8, -280 * y8, 0.43, 0.69, -1.81, 1.81, -280 * y6, 280 * y6, ...
11-
-280 * y6]);
12-
13-
% J2 = [
14-
% -1.71, 0.43, 8.32, 0, 0, 0, 0, 0;
15-
% 1.71, -8.75, 0, 0, 0, 0, 0, 0;
16-
% 0, 0, -10.03, 0.43, 0.035, 0, 0, 0;
17-
% 0, 8.32, 1.71, -1.12, 0, 0, 0, 0;
18-
% 0, 0, 0, 0, -1.745, 0.43, 0.43, 0;
19-
% 0, 0, 0, 0.69, 1.71, -280 * y8 - 0.43, 0.69, -280 * y6;
20-
% 0, 0, 0, 0, 0, 280 * y8, -1.81, 280 * y6;
21-
% 0, 0, 0, 0, 0, -280 * y8, 1.81, -280 * y6];
7+
jac = sparse( ...
8+
[1, 2, 1, 2, 4, 1, 3, 4, 3, 4, 6, 3, 5, 6, 5, 6, 7, 8, 5, 6, 7, 8, 6, 7, 8], ...
9+
[1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8], ...
10+
[-k1, k1, k2, -k2 - k3, k3, k6, -k1 - k6, k1, k2, -k2 - k4, k4, k5, -k1 - k5, k1, k2, -k2 - kPlusE, kPlusE, ...
11+
-kPlusE, k2, kMinus, -kSum, kSum, -kPlusPfrXPrime, kPlusPfrXPrime, -kPlusPfrXPrime]);
2212

2313
end

‎toolbox/+otp/+hires/jacobianAdjointVectorProduct.m

+13-21
Original file line numberDiff line numberDiff line change
@@ -1,25 +1,17 @@
1-
function Jv = jacobianAdjointVectorProduct(~, y, v)
1+
function vj = jacobianAdjointVectorProduct(~, y, v, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, ~)
22

3-
y6 = conj(y(6));
4-
y8 = conj(y(8));
3+
kPlusPfrXPrime = conj(kPlus * y(6));
4+
kPlusE = conj(kPlus * y(8));
5+
kSum = conj(k2 + kMinus + kStar);
56

6-
v1 = v(1);
7-
v2 = v(2);
8-
v3 = v(3);
9-
v4 = v(4);
10-
v5 = v(5);
11-
v6 = v(6);
12-
v7 = v(7);
13-
v8 = v(8);
14-
15-
Jv = [
16-
1.71 * (v2 - v1);
17-
0.43 * v1 - 8.75 * v2 + 8.32 * v4;
18-
8.32 * v1 - 10.03 * v3 + 1.71 * v4;
19-
0.43 * v3 - 1.12 * v4 + 0.69 * v6;
20-
0.035 * v3 - 1.745 * v5 + 1.71 * v6;
21-
0.43 * (v5 - v6) + 280 * y8 * (v7 - v6 - v8);
22-
0.43 * v5 + 0.69 * v6 + 1.81 * (v8 - v7);
23-
280 * y6 * (v7 - v6 - v8)];
7+
vj = [ ...
8+
conj(k1) * (v(2, :) - v(1, :));
9+
conj(k2) * v(1, :) - conj(k2 + k3) * v(2, :) + conj(k3) * v(4, :);
10+
conj(k6) * v(1, :) - conj(k1 + k6) * v(3, :) + conj(k1) * v(4, :);
11+
conj(k2) * v(3, :) - conj(k2 + k4) * v(4, :) + conj(k4) * v(6, :);
12+
conj(k5) * v(3, :) - conj(k1 + k5) * v(5, :) + conj(k1) * v(6, :);
13+
conj(k2) * v(5, :) - (conj(k2) + kPlusE) * v(6, :) + kPlusE * (v(7, :) - v(8, :));
14+
conj(k2) * v(5, :) + conj(kMinus) * v(6, :) + kSum * (v(8, :) - v(7, :));
15+
kPlusPfrXPrime * (v(7, :) - v(6, :) - v(8, :))];
2416

2517
end

0 commit comments

Comments
 (0)