Kuramoto Sivashinky 1D PDE

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

@@ -0,0 +1,45 @@
+classdef Canonical < otp.kuramotosivashinsky.KuramotoSivashinskyProblem
+    % The Kuramoto-Sivashinky equation is a chaotic problem.
+    %
+    % In this particular discretization, we are applying a spectral
+    % method, therefore the boundary conditions will be chosen to be as cyclic
+    % on the domain [0, L]. Note that this is different from another typical
+    % domain of [-L, L]. The larger the L, the more interesting the problem is
+    % but the more points are required to do a good discretization. The current
+    % canonical implementation with the size, L, and  is used in
+    %
+    % Kassam, Aly-Khan, and Lloyd N. Trefethen.
+    % "Fourth-order time-stepping for stiff PDEs."
+    % SIAM Journal on Scientific Computing 26, no. 4 (2005): 1214-1233.
+    %
+
+    methods
+        function obj = Canonical(varargin)
+
+            p = inputParser;
+            addParameter(p, 'Size', 128);
+            addParameter(p, 'L', 32*pi);
+            parse(p, varargin{:});
+
+            s = p.Results;
+
+            N = s.Size;
+            L = s.L;
+            params = otp.kuramotosivashinsky.KuramotoSivashinskyParameters;
+            params.L = L;
+
+            % exclude the left boundary point as it is identical to the
+            % right boundary point
+            x = linspace(L / N, L, N).';
+
+            u0 = cospi(2 * x / L) .* (1 + sinpi(2 * x / L));
+
+            u0hat = fft(u0);
+
+            tspan = [0, 150];
+
+            obj = obj@otp.kuramotosivashinsky.KuramotoSivashinskyProblem(tspan, u0hat, params);
+
+        end
+    end
+end

src/+otp/+kuramotosivashinsky/+presets/HairerWanner.m

@@ -0,0 +1,38 @@
+classdef HairerWanner < otp.kuramotosivashinsky.KuramotoSivashinskyProblem
+
+    methods
+        function obj = HairerWanner(varargin)
+
+            p = inputParser;
+            addParameter(p, 'Size', 1024);
+            addParameter(p, 'L', 80 * pi);
+
+            parse(p, varargin{:});
+
+            s = p.Results;
+
+            N = s.Size;
+            L = s.L;
+
+            params.L = L;
+
+            % exclude the left boundary point as it is identical to the
+            % right boundary point
+            x = linspace(L/N, L, N).';
+
+            eta1 = min(x/L, 0.1 - x/L);
+            eta2 = 20*(x/L - 0.2).*(0.3 - x/L);
+            eta3 = min(x/L - 0.6, 0.7 - x/L);
+            eta4 = min(x/L - 0.9, 1 - x/L);
+
+            u0 = 16*max(0, max([eta1, eta2, eta3, eta4], [], 2));
+
+            u0hat = fft(u0);
+
+            tspan = [0, 100];
+
+            obj = obj@otp.kuramotosivashinsky.KuramotoSivashinskyProblem(tspan, u0hat, params);
+
+        end
+    end
+end

src/+otp/+kuramotosivashinsky/KuramotoSivashinskyParameters.m

@@ -0,0 +1,7 @@
+classdef KuramotoSivashinskyParameters
+    % Parameters for the Kuramoto Sivashinsky problem
+    properties
+        % Represents the length-scale of the problem
+        L
+    end
+end

src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m

@@ -0,0 +1,32 @@
+classdef KuramotoSivashinskyProblem < otp.Problem
+
+    methods
+        function obj = KuramotoSivashinskyProblem(timeSpan, y0, parameters)
+            obj@otp.Problem('Kuramoto-Sivashinsky', [], timeSpan, y0, parameters);
+        end
+    end
+
+    methods
+        function soly = convert2grid(~, soly)
+            soly = abs(ifft(soly));
+        end
+    end
+
+    methods (Access = protected)
+        function onSettingsChanged(obj)
+            N = obj.NumVars;
+
+            k = 2 * pi * [0:(N/2 - 1), 0, (-N/2 + 1):-1].' / obj.Parameters.L;
+            ik = 1i * k;
+            k24 = k.^2 - k.^4;
+
+            obj.RHS = otp.RHS(@(t, u) otp.kuramotosivashinsky.f(t, u, ik, k24), ...
+                'Jacobian', ...
+                @(t, u) otp.kuramotosivashinsky.jac(t,u, ik, k24), ...
+                'JacobianVectorProduct', ...
+                @(t, u, v) otp.kuramotosivashinsky.jvp(t, u, v, ik, k24), ...
+                'JacobianAdjointVectorProduct', ...
+                @(t, u, v) otp.kuramotosivashinsky.javp(t, u, v, ik, k24));
+        end
+    end
+end

src/+otp/+kuramotosivashinsky/f.m

@@ -0,0 +1,5 @@
+function ut = f(~, u, ik, k24)
+
+ut = k24 .* u - (ik/2) .* fft(ifft(u).^2);
+
+end

src/+otp/+kuramotosivashinsky/jac.m

@@ -0,0 +1,6 @@
+function j = jac(~, u, ik, k24)
+
+circulantU = toeplitz(u, [u(1); flipud(u(2:end))]);
+j = diag(k24) - ik / length(u) .* circulantU;
+
+end

src/+otp/+kuramotosivashinsky/javp.m

@@ -0,0 +1,5 @@
+function jv = javp(~, u, v, ik, k24)
+
+jv = k24 .* v + fft(conj(ifft(u)) .* ifft(ik .* v));
+
+end

src/+otp/+kuramotosivashinsky/jvp.m

@@ -0,0 +1,5 @@
+function jv = jvp(~, u, v, ik, k24)
+
+jv = k24 .* v - ik .* fft(ifft(u) .* ifft(v));
+
+end