n-dimensional m-derivative generalized finite difference and an implementation of Korteweg de Vries using it

src/+otp/+kortewegdevries/+presets/AscherMcLachlan.m

@@ -0,0 +1,58 @@
+classdef AscherMcLachlan < otp.kortewegdevries.KortewegDeVriesProblem
+    methods
+        function obj = AscherMcLachlan
+            N = 200;
+
+            L = 1;
+
+            mesh = linspace(-L, L, N + 1);
+            mesh = mesh(1:(end-1));
+            meshBC = [];
+            order = 14;
+
+            u0 = cos(pi*mesh).';
+
+            BC = @(t, x) 0*x;
+
+            theta = 0;
+            nu = -(2/3)*(10^-3);
+            alpha = -3/8;
+            rho = -1/10;
+
+            hfun = @(x, y) hfunfull(x, y, L);
+
+            dist = @(x, y) abs(hfun(x, y));
+            r = 2*2/N;
+            weightfun = @(mesh, meshi) otp.utils.rbf.gaussian(dist(mesh, meshi)/r);
+
+            params = otp.kortewegdevries.KortewegDeVriesParameters;
+            params.GFDMOrder = order;
+            params.GFDMMesh = mesh;
+            params.GFDMMeshBC = meshBC;
+            params.GFDMHFun = hfun;
+            params.GFDMWeightFun = weightfun;
+            params.BCFun = BC;
+            params.Theta = theta;
+            params.Nu = nu;
+            params.Alpha = alpha;
+            params.Rho = rho;
+
+            timespan = [0, 5];
+
+            obj = [email protected](timespan, u0, params);
+
+        end
+    end
+end
+
+function h = hfunfull(x, y, L)
+
+hs = [x - y; 2*L + x - y; -2*L + x - y];
+
+[~, I] = min(abs(hs));
+
+J = 1:numel(x);
+
+h = hs(sub2ind(size(hs), I, J));
+
+end

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

@@ -0,0 +1,62 @@
+classdef Canonical < otp.kortewegdevries.KortewegDeVriesProblem
+    methods
+        function obj = Canonical
+            N = 250;
+
+            L = 10;
+
+            mesh = linspace(-L, L, N + 1);
+            mesh = mesh(1:(end-1));
+            meshBC = [];
+            order = 12;
+
+            u0 = 6*(sech(mesh).^2).';
+
+            BC = @(t, x) 0*x;
+
+            theta = 0;
+            nu = -1;
+            alpha = -3;
+            rho = 0;
+
+            hfun = @(x, y) hfunfull(x, y, L);
+
+            dist = @(x, y) abs(hfun(x, y));
+            r = 0.2;
+            weightfun = @(mesh, meshi) otp.utils.rbf.gaussian(dist(mesh, meshi)/r);
+
+            params = otp.kortewegdevries.KortewegDeVriesParameters;
+            params.GFDMOrder = order;
+            params.GFDMMesh = mesh;
+            params.GFDMMeshBC = meshBC;
+            params.GFDMHFun = hfun;
+            params.GFDMWeightFun = weightfun;
+            params.BCFun = BC;
+            params.Theta = theta;
+            params.Nu = nu;
+            params.Alpha = alpha;
+            params.Rho = rho;
+
+            timespan = [0, 2];
+
+            obj = [email protected](timespan, u0, params);
+
+        end
+
+
+    end
+
+
+end
+
+function h = hfunfull(x, y, L)
+
+hs = [x - y; 2*L + x - y; -2*L + x - y];
+
+[~, I] = min(abs(hs));
+
+J = 1:numel(x);
+
+h = hs(sub2ind(size(hs), I, J));
+
+end

src/+otp/+kortewegdevries/KortewegDeVriesParameters.m

@@ -0,0 +1,27 @@
+classdef KortewegDeVriesParameters
+
+properties
+    % 
+    GFDMOrder
+    % 
+    GFDMMesh
+    %
+    GFDMMeshBC
+    %
+    GFDMHFun
+    %
+    GFDMWeightFun
+    %
+    BCFun
+    %
+    Theta
+    %
+    Nu
+    %
+    Alpha
+    %
+    Rho
+end
+
+
+end

src/+otp/+kortewegdevries/KortewegDeVriesProblem.m

@@ -0,0 +1,74 @@
+classdef KortewegDeVriesProblem < otp.Problem
+
+    methods
+        function obj = KortewegDeVriesProblem(timeSpan, y0, parameters)
+
+            [email protected]('Korteweg de Vries', [], ...
+                timeSpan, y0, parameters);
+
+        end
+    end
+
+    properties (SetAccess = private)
+        DistanceFunction
+        Invariant
+    end
+
+    methods (Access = protected)
+        function onSettingsChanged(obj)
+
+            mesh = obj.Parameters.GFDMMesh;
+            meshBC = obj.Parameters.GFDMMeshBC;
+            order = obj.Parameters.GFDMOrder;
+            weightfun = obj.Parameters.GFDMWeightFun;
+            hfun = obj.Parameters.GFDMHFun;
+            BC = obj.Parameters.BCFun;
+            theta = obj.Parameters.Theta;
+            alpha = obj.Parameters.Alpha;
+            nu = obj.Parameters.Nu;
+            rho = obj.Parameters.Rho;
+
+            gridPts = size(mesh, 2);
+
+            if obj.NumVars ~= gridPts
+                warning('OTP:inconsistentNumVars', ...
+                    'NumVars is %d, but there are %d grid points', ...
+                    obj.NumVars, gridPts);
+            end
+
+            [A, B] = otp.utils.pde.gfdm.getAB(mesh, meshBC, order, weightfun, hfun);
+
+            obj.RHS = otp.RHS(@(t, u) ...
+                otp.kortewegdevries.f(t, u, BC, meshBC, A, B, theta, alpha, nu, rho));
+
+            %% Distance function and Invariant
+            dist = @(x, y) abs(hfun(x, y));
+
+            obj.DistanceFunction = @(t, y, i, j) dist(mesh(i), mesh(j));
+
+            % find the volume of the full mesh
+            meshfull = [mesh, meshBC];
+            volume = 0;
+            for i = 1:numel(mesh)
+                ds = sort(dist(meshfull, mesh(i)), "ascend");
+                volume = volume + sum(ds(1:3))/2;
+            end
+
+            obj.Invariant = @(t, u) ...
+                otp.kortewegdevries.invariant(t, u, BC, meshBC, A, B, theta, alpha, nu, rho, volume);
+
+        end
+
+        function sol = internalSolve(obj, varargin)
+            sol = [email protected](obj, 'Solver', otp.utils.Solver.Nonstiff, varargin{:});
+        end
+
+        function mov = internalMovie(obj, t, y, varargin)
+            % BUG This movie is technically broken, and will not work for
+            % arbitrary meshes.
+            mov = otp.utils.movie.LineMovie('title', obj.Name, varargin{:});
+            mov.record(t, y);
+        end
+
+    end
+end

src/+otp/+kortewegdevries/f.m

@@ -0,0 +1,16 @@
+function dudt = f(t, u, BC, meshBC, A, B, theta, alpha, nu, rho)
+
+uBC = BC(t, meshBC);
+
+dud = otp.utils.pde.gfdm.evalAB(A, B, u.', uBC);
+
+dudx = dud(1, :).';
+dudxxx = dud(3, :).';
+
+du2d = otp.utils.pde.gfdm.evalAB(A, B, (u.^2).', uBC);
+
+du2dx = du2d(1, :).';
+
+dudt = alpha*(theta*du2dx + 2*(1 - theta)*(u.*dudx)) + rho*dudx + nu*dudxxx;
+
+end

src/+otp/+kortewegdevries/invariant.m

@@ -0,0 +1,14 @@
+function c = invariant(t, u, BC, meshBC, A, B, ~, alpha, nu, rho, volume)
+
+c = zeros(3, size(u, 2));
+
+c(1, :) = volume*mean(u, 1);
+c(2, :) = volume*mean(u.^2, 1);
+
+uBC = BC(t, meshBC);
+dud = otp.utils.pde.gfdm.evalAB(A, B, u.', uBC);
+dudx = dud(1, :).';
+
+c(3, :) = volume*mean((alpha/3)*(u.^3) + (rho/2)*(u.^2) - (nu/2)*(dudx.^2), 1);
+
+end

src/+otp/+utils/+compatibility/pagemfun.m

@@ -0,0 +1,41 @@
+function X = pagemfun(A, transpA, B, transpB, mfun)
+
+if ndims(A) > 3 || ndims(B) > 3
+    error('OTP:notImplemented', 'Dimensions greater than 3 are not supported.');
+end
+
+% repmat A or B depending on the third dimension 
+N = max(size(A, 3), size(B, 3));
+
+if size(A, 3) == 1
+    repmat(A, 1, 1, N);
+elseif size(A, 3) ~= N
+    error('OTP:dimensionMismatch', 'The dimensions of A and B do not match.');
+end
+
+if size(B, 3) == 1
+    repmat(B, 1, 1, N);
+elseif size(B, 3) ~= N
+    error('OTP:dimensionMismatch', 'The dimensions of A and B do not match.');
+end
+
+% transpose A and B as required
+switch transpA 
+    case 'transpose'
+        A = permute(A, [2, 1, 3]);
+    case 'ctranspose'
+        A = permute(conj(A), [2, 1, 3]);
+end
+
+switch transpB
+    case 'transpose'
+        B = permute(B, [2, 1, 3]);
+    case 'ctranspose'
+        B = permute(conj(B), [2, 1, 3]);
+end
+
+for i = N:-1:1
+    X(:, :, i) = mfun(A(:, :, i), B(:, :, i));
+end
+
+end

src/+otp/+utils/+compatibility/pagemldivide.m

@@ -0,0 +1,21 @@
+function X = pagemldivide(A, B, C, D)
+% OCTAVE FIX: This function is a naive implementation of the pagemldivide
+% functionality for the purposes of supporting octave and legacy matlab
+% installations
+
+if nargin == 2
+    transpA = 'none';
+    transpB = 'none';
+elseif nargin == 4
+    transpA = B;
+    transpB = D;
+    B = C;
+else
+    error('OTP:invalidNumberOfArguments', 'The number of arguments has to be 2 or 4.');
+end
+
+mfun = @mldivide;
+
+X = otp.utils.compatibility.pagemfun(A, transpA, B, transpB, mfun);
+
+end

src/+otp/+utils/+compatibility/pagemtimes.m

@@ -14,42 +14,8 @@
     error('OTP:invalidNumberOfArguments', 'The number of arguments has to be 2 or 4.');
 end
 
-if ndims(A) > 3 || ndims(B) > 3
-    error('OTP:notImplemented', 'Dimensions greater than 3 are not supported.');
-end
-
-% repmat A or B depending on the third dimension 
-N = max(size(A, 3), size(B, 3));
-
-if size(A, 3) == 1
-    repmat(A, 1, 1, N);
-elseif size(A, 3) ~= N
-    error('OTP:dimensionMismatch', 'The dimensions of A and B do not match.');
-end
-
-if size(B, 3) == 1
-    repmat(B, 1, 1, N);
-elseif size(B, 3) ~= N
-    error('OTP:dimensionMismatch', 'The dimensions of A and B do not match.');
-end
+mfun = @mtimes;
 
-% transpose A and B as required
-switch transpA 
-    case 'transpose'
-        A = permute(A, [2, 1, 3]);
-    case 'ctranspose'
-        A = permute(conj(A), [2, 1, 3]);
-end
-
-switch transpB
-    case 'transpose'
-        B = permute(B, [2, 1, 3]);
-    case 'ctranspose'
-        B = permute(conj(B), [2, 1, 3]);
-end
-
-for i = N:-1:1
-    X(:, :, i) = A(:, :, i)*B(:, :, i);
-end
+X = otp.utils.compatibility.pagemfun(A, transpA, B, transpB, mfun);
 
 end

src/+otp/+utils/+pde/+gfdm/evalAB.m

@@ -0,0 +1,25 @@
+function dfd = evalAB(A, B, f, fBC)
+
+% OCTAVE FIX: find out if the functions pagemtimes and pagemldivide exist, and
+% if they does not, replace them with a compatible function
+if exist('pagemldivide', 'builtin') == 0
+    pmld = @otp.utils.compatibility.pagemldivide;
+else
+    pmld = @pagemldivide;
+end
+if exist('pagemtimes', 'builtin') == 0
+    pmt = @otp.utils.compatibility.pagemtimes;
+else
+    pmt = @pagemtimes;
+end
+
+ndivterms = size(A, 1);
+
+ffull = [f, fBC];
+
+fdiff = permute(-f + ffull.', [3, 1, 2]);
+Bfdiff = (pmt(B, 'none', fdiff, 'transpose'));
+
+dfd = reshape(pmld(A, Bfdiff), ndivterms, []);
+
+end