Start the Sanz-Serna Problem (#13)

* Start the Sanz-Serna Problem

* Fix PR comments

* Add linear-forcing splitting of RHS

* Use size(,1) instead of length

* Add partial derivative wrt time

* Make Jacobians matrices instead of function handles

Co-authored-by: Steven Roberts <sroberts994@gmail.com>

Author: reid-g

src/+otp/+sanzserna/+presets/Canonical.m

@@ -0,0 +1,12 @@
+classdef Canonical < otp.sanzserna.SanzSernaProblem
+    methods
+        function obj = Canonical(numGridCells)
+            if nargin < 1
+                numGridCells = 32;
+            end
+            
+            y0 = linspace(1/numGridCells + 1, 2, numGridCells).';
+            obj = obj@otp.sanzserna.SanzSernaProblem([0, 1], y0, struct);
+        end
+    end
+end

src/+otp/+sanzserna/SanzSernaProblem.m

@@ -0,0 +1,43 @@
+classdef SanzSernaProblem < otp.Problem
+    % This test problem comes from
+    % Sanz-Serna, J. M., Verwer, J. G., & Hundsdorfer, W. H. (1986). Convergence and order reduction of Runge-Kutta schemes
+    % applied to evolutionary problems in partial differential equations. Numerische Mathematik, 50(4), 405–418.
+    % https://doi.org/10.1007/BF01396661
+    
+    properties (SetAccess = private)
+        RhsLinear
+        RhsForcing
+    end
+    
+    methods
+        function obj = SanzSernaProblem(timeSpan, y0, parameters)
+            obj@otp.Problem('Sanz-Serna', [], timeSpan, y0, parameters);
+        end
+    end
+    
+    methods (Access=protected)
+        
+        function onSettingsChanged(obj)
+            
+            D = spdiags(ones(obj.NumVars, 1) * [obj.NumVars, -obj.NumVars], [-1, 0], obj.NumVars, obj.NumVars);
+            
+            x = linspace(1/obj.NumVars, 1, obj.NumVars).';
+            
+            obj.Rhs = otp.Rhs(@(t, y) otp.sanzserna.f(t, y, D, x), ...
+                otp.Rhs.FieldNames.Jacobian, otp.sanzserna.jac(D, x), ...
+                otp.Rhs.FieldNames.JacobianVectorProduct, @(t, y, v) otp.sanzserna.jvp(t, y, v, D, x), ...
+                otp.Rhs.FieldNames.PartialDerivativeTime, @(t, y) otp.sanzserna.pdt(t, y, D, x));
+            
+            obj.RhsLinear = otp.Rhs(@(t, y) otp.sanzserna.flinear(t, y, D, x), ...
+                otp.Rhs.FieldNames.Jacobian, otp.sanzserna.jaclinear(D, x));
+            obj.RhsForcing = otp.Rhs(@(t, y) otp.sanzserna.fforcing(t, y, D, x), ...
+                otp.Rhs.FieldNames.Jacobian, otp.sanzserna.jacforcing(D, x));
+        end
+        
+        function sol = internalSolve(obj, varargin)
+            sol = internalSolve@otp.Problem(obj, 'Method', @ode45, varargin{:});
+        end
+        
+        %true solution
+    end
+end

src/+otp/+sanzserna/f.m

@@ -0,0 +1,12 @@
+function dy = f(t, y, D, x)
+
+hInv = size(y, 1);
+
+forcingTerm = (t - x) / (1 + t)^2;
+
+forcingTerm(1) = forcingTerm(1) + (hInv / (1 + t));
+
+dy = D * y + forcingTerm;
+
+end
+

src/+otp/+sanzserna/fforcing.m

@@ -0,0 +1,10 @@
+function dy = fforcing(t, y, ~, x)
+
+hInv = size(y, 1);
+
+dy = (t - x) / (1 + t)^2;
+
+dy(1) = dy(1) + (hInv / (1 + t));
+
+end
+

src/+otp/+sanzserna/flinear.m

@@ -0,0 +1,6 @@
+function dy = flinear(~, y, D, ~)
+
+dy = D * y;
+
+end
+

src/+otp/+sanzserna/jac.m

@@ -0,0 +1,4 @@
+function D = jac(D, ~)
+
+end
+

src/+otp/+sanzserna/jacforcing.m

@@ -0,0 +1,7 @@
+function J = jacforcing(~, x)
+
+numVars = length(x);
+J = sparse(numVars, numVars);
+
+end
+

src/+otp/+sanzserna/jaclinear.m

@@ -0,0 +1,4 @@
+function D = jaclinear(D, ~)
+
+end
+

src/+otp/+sanzserna/jvp.m

@@ -0,0 +1,6 @@
+function jx = jvp(~, ~, v, D, ~)
+
+jx = D * v;
+
+end
+

src/+otp/+sanzserna/pdt.m

@@ -0,0 +1,9 @@
+function jt = pdt(t, y, ~, x)
+
+hInv = size(y, 1);
+
+jt = (1 - t + 2 * x) / (1 + t)^3;
+jt(1) = jt(1) - hInv / (1 + t)^2;
+
+end
+