ComputationalScienceLaboratory · Steven-Roberts · Jan 28, 2025 · Oct 13, 2024 · Oct 14, 2024
diff --git a/docs/references.bib b/docs/references.bib
@@ -173,9 +173,39 @@ @article{Pet86
  number = {4},
  pages = {837--852},
  publisher = {Society for Industrial and Applied Mathematics},
- title = {Order Results for Implicit Runge-Kutta Methods Applied to Differential/ Algebraic Systems},
+ title = {Order Results for Implicit {Runge--Kutta} Methods Applied to Differential/Algebraic Systems},
  volume = {23},
  year = {1986},
- url = {https://www.jstor.org/stable/2157625}
+ doi = {10.1137/0723054}
+}
+
+@article{Sch75,
+  title = {A new approach to explain the ``high irradiance responses'' of photomorphogenesis on the basis of phytochrome},
+  volume = {2},
+  ISSN = {1432-1416},
+  DOI = {10.1007/bf00276015},
+  number = {1},
+  journal = {Journal of Mathematical Biology},
+  author = {Schäfer,  E.},
+  year = {1975},
+  pages = {41-56}
+}
+
+@article{Got77,
+  title = {{MISS} - ein einfaches Simulations-System für biologische und chemische Prozesse},
+  volume = {3},
+  journal = {EDV in Medizin und Biologie},
+  author = {B. A. Gottwald},
+  year = {1977},
+  pages = {85-90}
+}
+
+@article{dSL98,
+  title={Collecting real-life problems to test solvers for implicit differential equations},
+  author={de Swart, Jacques J. B. and Lioen, Walter M.},
+  journal={CWI Quarterly},
+  volume={11},
+  number={1},
+  pages={83-100},
+  year={1998}
 }
-
diff --git a/toolbox/+otp/+hires/+presets/Canonical.m b/toolbox/+otp/+hires/+presets/Canonical.m
@@ -1,9 +1,41 @@
 classdef Canonical < otp.hires.HIRESProblem
+    % HIRES preset proposed in :cite:p:`HW96` (pp. 144-145) which uses time span $t ∈ [0, 321.8122]$, initial conditions
+    %
+    % $$
+    % y(0) = [1, 0, 0, 0, 0, 0, 0, 0.0057]^T,
+    % $$
+    %
+    % and parameters
+    %
+    % $$
+    % k_1 &= 1.71, \quad & k_2 &= 0.43, \quad & k_3 &= 8.32, \quad & k_4 &= 0.69, \quad & k_5 &= 0.035, \\
+    % k_6 &= 8.32, & k_{+} &= 280, & k_{-} &= 0.69, & k^* &= 0.69, & o_{k_s} &= 0.0007.
+    % $$
+
     methods
         function obj = Canonical(varargin)
+            % Create the canonical HIRES problem object.
+            %
+            % Parameters
+            % ----------
+            % varargin
+            %    A variable number of name-value pairs. The accepted names are
+            %
+            %    - ``K1`` – Value of $k_1$.
+            %    - ``K2`` – Value of $k_2$.
+            %    - ``K3`` – Value of $k_3$.
+            %    - ``K4`` – Value of $k_4$.
+            %    - ``K5`` – Value of $k_5$.
+            %    - ``K6`` – Value of $k_6$.
+            %    - ``KPlus`` – Value of $k_{+}$.
+            %    - ``KMinus`` – Value of $k_{-}$.
+            %    - ``KStar`` – Value of $k^*$.
+            %    - ``OKS`` – Value of $o_{k_s}$.
+
             tspan = [0, 321.8122];
             y0 = [1; 0; 0; 0; 0; 0; 0; 0.0057];
-            params = otp.hires.HIRESParameters(varargin{:});
+            params = otp.hires.HIRESParameters('K1', 1.71, 'K2', 0.43, 'K3', 8.32, 'K4', 0.69, 'K5', 0.035, ...
+                'K6', 8.32, 'KPlus', 280, 'KMinus', 0.69, 'KStar', 0.69, 'OKS', 0.0007, varargin{:});
             obj = obj@otp.hires.HIRESProblem(tspan, y0, params);
         end
     end

diff --git a/toolbox/+otp/+hires/HIRESParameters.m b/toolbox/+otp/+hires/HIRESParameters.m
@@ -1,6 +1,38 @@
 classdef HIRESParameters < otp.Parameters
     % Parameters for the HIRES problem.
 
+    properties
+        % The reaction rate $k_1$.
+        K1 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical}
+
+        % The reaction rate $k_2$.
+        K2 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical}
+
+        % The reaction rate $k_3$.
+        K3 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical}
+
+        % The reaction rate $k_4$.
+        K4 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical}
+
+        % The reaction rate $k_5$.
+        K5 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical}
+
+        % The reaction rate $k_6$.
+        K6 %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical}
+
+        % The reaction rate $k_{+}$.
+        KPlus %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical}
+
+        % The reaction rate $k_{-}$.
+        KMinus %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical}
+
+        % The reaction rate $k^*$.
+        KStar %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical}
+
+        % The source term $o_{k_s}$.
+        OKS %MATLAB ONLY: (1,1) {otp.utils.validation.mustBeNumerical}
+    end
+
     methods
         function obj = HIRESParameters(varargin)
             % Create a HIRES parameters object.

diff --git a/toolbox/+otp/+hires/HIRESProblem.m b/toolbox/+otp/+hires/HIRESProblem.m
@@ -1,18 +1,117 @@
-
 classdef HIRESProblem < otp.Problem
+    % An eight-variable high irradiance response model from plant physiology.
+    %
+    % In :cite:p:`Sch75`, the following chemical reactions were proposed to explain the high irradiance responses of
+    % photomorphogenesis on the basis of phytochrome:
+    %
+    % $$
+    % \begin{array}{cccc}
+    % \xrightarrow{o_{k_s}} & \ce{P_{{r}}} & \xrightleftharpoons[k_2]{k_1} & \ce{P_{fr}} \\
+    % & {\scriptsize k_6} ↑ & & ↓ {\scriptsize k_3} \\
+    % & \ce{P_{{r}}X} & \xrightleftharpoons[k_2]{k_1} & \ce{P_{fr}X} \\
+    % & {\scriptsize k_5} ↑ & & ↓ {\scriptsize k_4} \\
+    % & \ce{P_{{r}}X'} & \xrightleftharpoons[k_2]{k_1} & \ce{P_{fr}X'} \\
+    % \end{array}
+    % $$
+    %
+    % $$
+    % \begin{array}{ccccc}
+    % \ce{E + P_{{r}}X'} & \xleftarrow{k_2} & \ce{P_{fr}X'E} & \xrightleftharpoons[k_{-}]{k_{+}} & \ce{P_{fr}X' + E} \\
+    % & & ↓ {\scriptsize k^*} \\
+    % & & \ce{P_{fr}' + X' + E}
+    % \end{array}
+    % $$
+    %
+    % Reactants $\ce{P_{{r}}}$ and $\ce{P_{fr}}$ are the red and far-red absorbing form of phytochrome, respectively.
+    % These can be bound by receptors $\ce{X}$ and $\ce{X'}$, partially influenced by enzyme $\ce{E}$. The system is
+    % modeled by the differential equations :cite:p:`Got77` :cite:p:`dSL98`
+    %
+    % $$
+    % \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}, \\
+    % \frac{d}{dt} \ce{[P_{fr}]} &= k_1 \ce{[P_{{r}}]} - (k_2 + k_3) \ce{[P_{fr}]}, \\
+    % \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']}, \\
+    % \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]}, \\
+    % \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]}, \\
+    % \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']}
+    % + k_{-} \ce{[P_{fr}X'E]} - k_{+} \ce{[P_{fr}X']} \ce{[E]}, \\
+    % \frac{d}{dt} \ce{[P_{fr}X'E]} &= -(k_2 + k_{-} + k^*) \ce{[P_{fr}X'E]} + k_{+} \ce{[P_{fr}X']} \ce{[E]}, \\
+    % \frac{d}{dt} \ce{[E]} &= (k_2 + k_{-} + k^*) \ce{[P_{fr}X'E]} - k_{+} \ce{[P_{fr}X']} \ce{[E]}.
+    % $$
+    %
+    % Notes
+    % -----
+    % +---------------------+--------------------------------------------------------------------+
+    % | Type                | ODE                                                                |
+    % +---------------------+--------------------------------------------------------------------+
+    % | Number of Variables | 8                                                                  |
+    % +---------------------+--------------------------------------------------------------------+
+    % | Stiff               | typically, depending on $k_1, …, k_6$, $k_{+}$, $k_{-}$, and $k^*$ |
+    % +---------------------+--------------------------------------------------------------------+
+    %
+    % Example
+    % -------
+    % >>> problem = otp.hires.presets.Canonical;
+    % >>> sol = problem.solve('AbsTol', 1e-12);
+    % >>> problem.plot(sol);
 
     methods
         function obj = HIRESProblem(timeSpan, y0, parameters)
+            % Create a HIRES problem object.
+            %
+            % Parameters
+            % ----------
+            % timeSpan : numeric(1, 2)
+            %    The start and final time.
+            % y0 : numeric(2, 1)
+            %    The initial conditions.
+            % parameters : otp.hires.HIRESParameters
+            %    The parameters.
             obj@otp.Problem('HIRES', 8, timeSpan, y0, parameters);
         end
     end
 
     methods (Access = protected)
         function onSettingsChanged(obj)
-            obj.RHS = otp.RHS(@otp.hires.f, ...
-                'Jacobian', @otp.hires.jacobian, ...
-                'JacobianVectorProduct', @otp.hires.jacobianVectorProduct, ...
-                'JacobianAdjointVectorProduct', @otp.hires.jacobianAdjointVectorProduct);
+            k1 = obj.Parameters.K1;
+            k2 = obj.Parameters.K2;
+            k3 = obj.Parameters.K3;
+            k4 = obj.Parameters.K4;
+            k5 = obj.Parameters.K5;
+            k6 = obj.Parameters.K6;
+            kPlus = obj.Parameters.KPlus;
+            kMinus = obj.Parameters.KMinus;
+            kStar = obj.Parameters.KStar;
+            oks = obj.Parameters.OKS;
+
+            obj.RHS = otp.RHS(@(t, y) otp.hires.f(t, y, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, oks), ...
+                'Jacobian', @(t, y) otp.hires.jacobian(t, y, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, oks), ...
+                'JacobianVectorProduct', @(t, y, v) ...
+                otp.hires.jacobianVectorProduct(t, y, v, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, oks), ...
+                'JacobianAdjointVectorProduct', @(t, y, v) ...
+                otp.hires.jacobianAdjointVectorProduct(t, y, v, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, oks), ...
+                'Vectorized', 'on', ...
+                'NonNegative', 1:obj.NumVars);
+        end
+
+        function label = internalIndex2label(obj, index)
+            switch  index
+                case 1
+                    label = 'P_r';
+                case 2
+                    label = 'P_{fr}';
+                case 3
+                    label = 'P_r X';
+                case 4
+                    label = 'P_{fr} X';
+                case 5
+                    label = 'P_r X''';
+                case 6
+                    label = 'P_{fr} X''';
+                case 7
+                    label = 'P_{fr} X'' E';
+                case 8
+                    label = 'E';
+            end
         end
 
         function fig = internalPlot(obj, t, y, varargin)

diff --git a/toolbox/+otp/+hires/f.m b/toolbox/+otp/+hires/f.m
@@ -1,13 +1,28 @@
-function yPrime = f(~, y)
+function yPrime = f(~, y, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, oks)
 
-yPrime = [
-    -1.71 * y(1) + 0.43 * y(2) + 8.32 * y(3) + 0.0007;
-    1.71 * y(1) - 8.75 * y(2);
-    -10.03 * y(3) + 0.43 * y(4) + 0.035 * y(5);
-    8.32 * y(2) + 1.71 * y(3) - 1.12 * y(4);
-    -1.745 * y(5) + 0.43 * y(6) + 0.43 * y(7);
-    -280 * y(6) * y(8) + 0.69 * y(4) + 1.71 * y(5) - 0.43 * y(6) + 0.69 * y(7);
-    280 * y(6) * y(8) - 1.81 * y(7);
-    -280 * y(6) * y(8) + 1.81 * y(7)];
+k1Pr = k1 * y(1, :);
+k2Pfr = k2 * y(2, :);
+k3Pfr = k3 * y(2, :);
+k1PrX = k1 * y(3, :);
+k6PrX = k6 * y(3, :);
+k2PfrX = k2 * y(4, :);
+k4PfrX = k4 * y(4, :);
+k1PrXPrime = k1 * y(5, :);
+k5PrXPrime = k5 * y(5, :);
+k2PfrXPrime = k2 * y(6, :);
+k2PfrXPrimeE = k2 * y(7, :);
+kMinusPfrXPrimeE = kMinus * y(7, :);
+kPlusPfrXPrimeE = kPlus * y(6, :) * y(8, :);
+yPrime7 = -k2PfrXPrimeE - (kMinus + kStar) * y(7, :) + kPlusPfrXPrimeE;
 
-end
+yPrime = [ ...
+    -k1Pr + k2Pfr + k6PrX + oks;
+    k1Pr - k2Pfr - k3Pfr;
+    -k1PrX - k6PrX + k2PfrX + k5PrXPrime;
+    k3Pfr + k1PrX - k2PfrX - k4PfrX;
+    -k1PrXPrime - k5PrXPrime + k2PfrXPrime + k2PfrXPrimeE;
+    k4PfrX + k1PrXPrime - k2PfrXPrime + kMinusPfrXPrimeE - kPlusPfrXPrimeE;
+    yPrime7;
+    -yPrime7];
+
+end
diff --git a/toolbox/+otp/+hires/jacobian.m b/toolbox/+otp/+hires/jacobian.m
@@ -1,23 +1,13 @@
-function J = jacobian(~, y)
+function jac = jacobian(~, y, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, ~)
 
-y6 = y(6);
-y8 = y(8);
+kPlusPfrXPrime = kPlus * y(6);
+kPlusE = kPlus * y(8);
+kSum = k2 + kMinus + kStar;
 
-J = sparse( ...
-    [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], ...
-    [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], ...
-    [-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, ...
-    0.43, -280 * y8 - 0.43, 280 * y8, -280 * y8, 0.43, 0.69, -1.81, 1.81, -280 * y6, 280 * y6, ...
-    -280 * y6]);
-
-% J2 = [
-%     -1.71, 0.43, 8.32, 0, 0, 0, 0, 0;
-%     1.71, -8.75, 0, 0, 0, 0, 0, 0;
-%     0, 0, -10.03, 0.43, 0.035, 0, 0, 0;
-%     0, 8.32, 1.71, -1.12, 0, 0, 0, 0;
-%     0, 0, 0, 0, -1.745, 0.43, 0.43, 0;
-%     0, 0, 0, 0.69, 1.71, -280 * y8 - 0.43, 0.69, -280 * y6;
-%     0, 0, 0, 0, 0, 280 * y8, -1.81, 280 * y6;
-%     0, 0, 0, 0, 0, -280 * y8, 1.81, -280 * y6];
+jac = sparse( ...
+    [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], ...
+    [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], ...
+    [-k1, k1, k2, -k2 - k3, k3, k6, -k1 - k6, k1, k2, -k2 - k4, k4, k5, -k1 - k5, k1, k2, -k2 - kPlusE, kPlusE, ...
+    -kPlusE, k2, kMinus, -kSum, kSum, -kPlusPfrXPrime, kPlusPfrXPrime, -kPlusPfrXPrime]);
 
 end
diff --git a/toolbox/+otp/+hires/jacobianAdjointVectorProduct.m b/toolbox/+otp/+hires/jacobianAdjointVectorProduct.m
@@ -1,25 +1,17 @@
-function Jv = jacobianAdjointVectorProduct(~, y, v)
+function vj = jacobianAdjointVectorProduct(~, y, v, k1, k2, k3, k4, k5, k6, kPlus, kMinus, kStar, ~)
 
-y6 = conj(y(6));
-y8 = conj(y(8));
+kPlusPfrXPrime = conj(kPlus * y(6));
+kPlusE = conj(kPlus * y(8));
+kSum = conj(k2 + kMinus + kStar);
 
-v1 = v(1);
-v2 = v(2);
-v3 = v(3);
-v4 = v(4);
-v5 = v(5);
-v6 = v(6);
-v7 = v(7);
-v8 = v(8);
-
-Jv = [
-    1.71 * (v2 - v1);
-    0.43 * v1 - 8.75 * v2 + 8.32 * v4;
-    8.32 * v1 - 10.03 * v3 + 1.71 * v4;
-    0.43 * v3 - 1.12 * v4 + 0.69 * v6;
-    0.035 * v3 - 1.745 * v5 + 1.71 * v6;
-    0.43 * (v5 - v6) + 280 * y8 * (v7 - v6 - v8);
-    0.43 * v5 + 0.69 * v6 + 1.81 * (v8 - v7);
-    280 * y6 * (v7 - v6 - v8)];
+vj = [ ...
+    conj(k1) * (v(2, :) - v(1, :));
+    conj(k2) * v(1, :) - conj(k2 + k3) * v(2, :) + conj(k3) * v(4, :);
+    conj(k6) * v(1, :) - conj(k1 + k6) * v(3, :) + conj(k1) * v(4, :);
+    conj(k2) * v(3, :) - conj(k2 + k4) * v(4, :) + conj(k4) * v(6, :);
+    conj(k5) * v(3, :) - conj(k1 + k5) * v(5, :) + conj(k1) * v(6, :);
+    conj(k2) * v(5, :) - (conj(k2) + kPlusE) * v(6, :) + kPlusE * (v(7, :) - v(8, :));
+    conj(k2) * v(5, :) + conj(kMinus) * v(6, :) + kSum * (v(8, :) - v(7, :));
+    kPlusPfrXPrime * (v(7, :) - v(6, :) - v(8, :))];
 
 end