diff --git a/src/+otp/+rangangermann/+presets/Canonical.m b/src/+otp/+rangangermann/+presets/Canonical.m new file mode 100644 index 00000000..86c1176c --- /dev/null +++ b/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 diff --git a/src/+otp/+rangangermann/RangAngermannParameters.m b/src/+otp/+rangangermann/RangAngermannParameters.m new file mode 100644 index 00000000..bdc68cfb --- /dev/null +++ b/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 diff --git a/src/+otp/+rangangermann/RangAngermannProblem.m b/src/+otp/+rangangermann/RangAngermannProblem.m new file mode 100644 index 00000000..e7022a34 --- /dev/null +++ b/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 diff --git a/src/+otp/+rangangermann/f.m b/src/+otp/+rangangermann/f.m new file mode 100644 index 00000000..9ad47a25 --- /dev/null +++ b/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 diff --git a/src/+otp/+rangangermann/jacobian.m b/src/+otp/+rangangermann/jacobian.m new file mode 100644 index 00000000..829e85d3 --- /dev/null +++ b/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 diff --git a/src/+otp/+rangangermann/mass.m b/src/+otp/+rangangermann/mass.m new file mode 100644 index 00000000..95452e94 --- /dev/null +++ b/src/+otp/+rangangermann/mass.m @@ -0,0 +1,5 @@ +function M = mass(n) + + + +end