Rang Angermann problem

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

@@ -0,0 +1,43 @@
+classdef Canonical < otp.rangangermann.RangAngermannProblem
+    %CANONICAL
+
+    methods
+        function obj = Canonical(varargin)
+
+            p = inputParser;
+            p.addParameter('Size', 256);
+
+            % exact solution
+            exactu = @(t, x, y) (2*x + y).*sin(t);
+            exactv = @(t, x, y) (x + 3*y).*cos(t);
+
+            p.parse(varargin{:});
+
+            s = p.Results;
+            n = s.Size;
+
+            params = otp.rangangermann.RangAngermannParameters;
+            params.Size = n;
+            % the forcing functions are based on the excact solution
+            params.ForcingU = @(t, x, y) (2*x + y).*cos(t) + 2*x.*sin(t) + y.*sin(t) ...
+                - exactu(t, x, y) + exactv(t, x, y);
+            params.ForcingV = @(t, x, y) exactu(t, x, y).^2 + exactv(t, x, y).^2;
+
+            % the boundary conditions are the exact solutions
+            params.BoundaryConditionsU = exactu;
+            params.BoundaryConditionsV = exactv;
+
+            x = linspace(0, 1, n + 2);
+            y = linspace(0, 1, n + 2);
+            [yfull, xfull] = meshgrid(x, y);
+
+            x = xfull(2:(end-1), 2:(end-1));
+            y = yfull(2:(end-1), 2:(end-1));
+            uv0 = [reshape(exactu(0, x, y), [], 1); reshape(exactv(0, x, y), [], 1)];
+
+            tspan = [0, 1];
+
+            obj = obj@otp.rangangermann.RangAngermannProblem(tspan, uv0, params);
+        end
+    end
+end

src/+otp/+rangangermann/RangAngermannParameters.m

@@ -0,0 +1,19 @@
+classdef RangAngermannParameters
+    % User-configurable parameters for the Rang-Angermann Problem
+    %
+    %   See also otp.rangangermann.RangAngermannProblem
+
+    properties
+        %Size is the dimension of the problem
+        Size %MATLAB ONLY: (1,1) {mustBeInteger}
+        %ForcingU is a forcing function
+        ForcingU %MATLAB ONLY: {mustBeA(ForcingU, {'numeric', 'function_handle'})}
+        %ForcingV is a forcing function
+        ForcingV %MATLAB ONLY: {mustBeA(ForcingU, {'numeric', 'function_handle'})}
+        %BoundaryConditionsU is a boundary
+        BoundaryConditionsU %MATLAB ONLY: {mustBeA(BoundaryConditionsU, {'numeric', 'function_handle'})}
+        %BoundaryConditionsV is a boundary
+        BoundaryConditionsV %MATLAB ONLY: {mustBeA(BoundaryConditionsV, {'numeric', 'function_handle'})}
+
+    end
+end

src/+otp/+rangangermann/RangAngermannProblem.m

@@ -0,0 +1,106 @@
+classdef RangAngermannProblem < otp.Problem
+    % This is an Index-1 PDAE
+    %
+    % See
+    %  Rang, J., & Angermann, L. (2005). 
+    %  New Rosenbrock W-methods of order 3 for partial differential algebraic equations of index 1. 
+    %  BIT Numerical Mathematics, 45, 761-787.
+    %
+
+    methods
+        function obj = RangAngermannProblem(timeSpan, y0, parameters)
+            obj@otp.Problem('Rang-Angermann Problem', [], timeSpan, y0, parameters);
+        end
+    end
+
+    methods (Access = protected)
+        function onSettingsChanged(obj)
+            n = obj.Parameters.Size;
+            forcingu = obj.Parameters.ForcingU;
+            forcingv = obj.Parameters.ForcingV;
+            FBCu = obj.Parameters.BoundaryConditionsU;
+            FBCv = obj.Parameters.BoundaryConditionsV;
+
+            if obj.NumVars ~= n^2
+                warning('OTP:inconsistentNumVars', ...
+                    'NumVars is %d, but there are %d grid points', ...
+                    obj.NumVars, n^2);
+            end
+
+            % First we construct the finite difference operators in the
+            % full domain including the boundary conditions
+            nfull = n + 2;
+
+            hx = 1/(n + 3);
+            hy = 1/(n + 3);
+            % one dimensional derivative in the x direciton
+            Ddx = spdiags(repmat([-1 1 0]/(hx), nfull, 1), [-1, 0, 1], nfull, nfull);
+            % two dimensional derivative in the x direction
+            Dx = kron(speye(nfull), Ddx);
+
+            % one dimensional derivative in the y direciton
+            Ddy = spdiags(repmat([-1 1 0]/(hy), nfull, 1), [-1, 0, 1], nfull, nfull);
+            % two dimensional derivative in the y direction
+            Dy = kron(Ddy, speye(nfull));
+
+            % both one dimensional Laplacians
+            Ldx = spdiags(repmat([1 -2 1]/(hx^2), nfull, 1), [-1, 0, 1], nfull, nfull);
+            Ldy = spdiags(repmat([1 -2 1]/(hy^2), nfull, 1), [-1, 0, 1], nfull, nfull);
+
+            % two dimensional Laplacian
+            L = kron(speye(nfull), Ldx) + kron(Ldy, speye(nfull));
+
+            % construct the x-y grid
+            xb = linspace(0, 1, nfull);
+            yb = linspace(0, 1, nfull);
+
+            [y, x] = meshgrid(xb, yb);
+
+
+            % internal point grid
+            xinternal = x(2:(end - 1), 2:(end - 1));
+            yinternal = y(2:(end - 1), 2:(end - 1));
+
+            % construct the x derivative operator on the internal grid
+            Dxinternal = kron(speye(n), Ddx(2:(end - 1), 2:(end - 1)));
+            % construct the y derivative operator on the internal grid
+            Dyinternal = kron(Ddy(2:(end - 1), 2:(end - 1)), speye(n));
+            % construct the Laplacian on the internal grid
+            Linternal = kron(speye(n), Ldx(2:(end - 1), 2:(end - 1))) ...
+                + kron(Ldy(2:(end - 1), 2:(end - 1)), speye(n));
+
+            % construct the mass matrix
+            Md = [ones(n^2, 1); zeros(n^2, 1)];
+            M = spdiags(Md, 0, 2*(n^2), 2*(n^2));
+
+            % construct the RHS
+            obj.RHS = otp.RHS(@(t, uv) otp.rangangermann.f(t, uv, L, Dx, Dy, x, y, forcingu, forcingv, FBCu, FBCv), ...
+                'Jacobian', @(t, uv) otp.rangangermann.jacobian(t, uv, Linternal, Dxinternal, Dyinternal, xinternal, yinternal, forcingu, forcingv, FBCu, FBCv), ...
+                'Mass', M);
+
+        end
+
+        function uv = internalSolveExactly(obj, t)
+            n = obj.Parameters.Size;
+            nfull = n + 2;
+
+            % construct the x-y grid
+            xb = linspace(0, 1, nfull);
+            yb = linspace(0, 1, nfull);
+
+            [y, x] = meshgrid(xb, yb);
+
+            % internal point grid
+            xinternal = x(2:(end - 1), 2:(end - 1));
+            yinternal = y(2:(end - 1), 2:(end - 1));
+
+            for i = length(t):-1:1
+                u = obj.Parameters.BoundaryConditionsU(t(i), xinternal, yinternal);
+                u = reshape(u, [], 1);
+                v = obj.Parameters.BoundaryConditionsV(t(i), xinternal, yinternal);
+                v = reshape(v, [], 1);
+                uv(:, i) = [u; v];
+            end
+        end
+    end
+end

src/+otp/+rangangermann/f.m

@@ -0,0 +1,39 @@
+function duvdt = f(t, uv, L, Dx, Dy, x, y, forcingu, forcingv, FBCu, FBCv)
+
+usmall = uv(1:(end/2));
+vsmall = uv((end/2 + 1):end);
+
+n = sqrt(numel(usmall));
+
+u = FBCu(t, x, y);
+v = FBCv(t, x, y);
+
+u(2:(end - 1), 2:(end - 1)) = reshape(usmall, n, n);
+v(2:(end - 1), 2:(end - 1)) = reshape(vsmall, n, n);
+
+u = reshape(u, [], 1);
+v = reshape(v, [], 1);
+
+Lu = L*u;
+Lv = L*v;
+Dxu = Dx*u;
+Dyu = Dy*u;
+
+xDxu = reshape(x, [], 1).*Dxu;
+yDyu = reshape(y, [], 1).*Dyu;
+
+f1 = reshape(forcingu(t, x, y), [], 1);
+f2 = reshape(forcingv(t, x, y), [], 1);
+
+dudt = f1 + Lu + Lv - xDxu - yDyu + u - v;
+valg = f2 + Lu + Lv - u.^2 - v.^2;
+
+dudt = reshape(dudt, n + 2, n + 2);
+dudt = reshape(dudt(2:(end - 1), 2:(end - 1)), [], 1);
+
+valg = reshape(valg, n + 2, n + 2);
+valg = reshape(valg(2:(end - 1), 2:(end - 1)), [], 1);
+
+duvdt = [dudt; valg];
+
+end

src/+otp/+rangangermann/jacobian.m

@@ -0,0 +1,19 @@
+function J = jacobian(~, uv, L, Dx, Dy, x, y, ~, ~, ~, ~)
+
+usmall = uv(1:(end/2));
+vsmall = uv((end/2 + 1):end);
+
+n = sqrt(numel(usmall));
+
+ddudu = L ...
+    - spdiags(reshape(x, [], 1), 0, n^2, n^2)*Dx ...
+    - spdiags(reshape(y, [], 1), 0, n^2, n^2)*Dy ...
+    + speye(n^2, n^2);
+ddudv = L - speye(n^2, n^2);
+ddvdu = L - spdiags(reshape(2*usmall, [], 1), 0, n^2, n^2);
+ddvdv = L - spdiags(reshape(2*vsmall, [], 1), 0, n^2, n^2);
+
+
+J = [ddudu, ddudv; ddvdu, ddvdv];
+
+end

src/+otp/+rangangermann/mass.m

@@ -0,0 +1,5 @@
+function M = mass(n)
+
+
+
+end