Kuramoto Sivashinky 1D PDE

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

@@ -0,0 +1,47 @@
+% 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.
+%
+
+classdef Canonical < otp.kuramotosivashinsky.KuramotoSivashinskyProblem
+
+    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.L = L;
+
+            h = L/N;
+
+            x = h*(1:N).';
+
+            u0 = cos(x/16).*(1+sin(x/16));
+
+            u0hat = fft(u0);
+
+            tspan = [0, 150];
+
+            obj = [email protected](tspan, u0hat, params);
+
+        end
+    end
+end

src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m

@@ -0,0 +1,57 @@
+classdef KuramotoSivashinskyProblem < otp.Problem
+
+    methods
+        function obj = KuramotoSivashinskyProblem(timeSpan, y0, parameters)
+            [email protected]('Kuramoto-Sivashinsky', [], timeSpan, y0, parameters);
+        end
+    end
+
+    methods
+        function soly = solution2real(soly)
+
+            if isstruct(soly)
+                soly.y = real(ifft(soly.y));
+            else
+                soly = real(ifft(soly.')).';
+            end
+
+        end
+
+    end
+
+    methods (Access = protected)
+        function onSettingsChanged(obj)
+
+            L = obj.Parameters.L;
+
+            N = obj.NumVars;
+
+            div = L/(2*pi);
+
+            k = (1i*[0:(N/2 - 1), 0, (-N/2 + 1):-1].'/div);
+            k2 = k.^2;
+            k4 = k.^4;
+
+            obj.Rhs = otp.Rhs(@(t, u) otp.kuramotosivashinsky.f(t, u, k, k2, k4), ...
+                otp.Rhs.FieldNames.Jacobian, ...
+                @(t, u) otp.kuramotosivashinsky.jac(t,u, k, k2, k4), ...
+                otp.Rhs.FieldNames.JacobianVectorProduct, ...
+                @(t, u, v) otp.kuramotosivashinsky.jvp(t,u, k, k2, k4, v), ...
+                otp.Rhs.FieldNames.JacobianAdjointVectorProduct, ...
+                @(t, u, v) otp.kuramotosivashinsky.javp(t,u, k, k2, k4, v));
+
+        end
+
+        function validateNewState(obj, newTimeSpan, newY0, newParameters)
+
+            [email protected](obj, ...
+                newTimeSpan, newY0, newParameters)
+
+            assert(mod(numel(newY0), 2) == 0);
+
+            otp.utils.StructParser(newParameters) ...
+                .checkField('L', 'scalar', 'finite', 'positive');
+
+        end
+    end
+end

src/+otp/+kuramotosivashinsky/f.m

@@ -0,0 +1,7 @@
+function ut = f(~, u, k, k2, k4)
+
+u2 = fft(real(ifft(u)).^2);
+
+ut = - k2.*u - k4.*u -(k/2).*u2;
+
+end

src/+otp/+kuramotosivashinsky/jac.m

@@ -0,0 +1,5 @@
+function j = jac(~, u, k, k2, k4)
+
+j = -diag(k2) - diag(k4) - diag(k)*ifft(fft(diag(ifft(u))).').';
+
+end

src/+otp/+kuramotosivashinsky/javp.m

@@ -0,0 +1,5 @@
+function jv = javp(~, u, k, k2, k4, v)
+
+jv = -k2.*v - k4.*v - fft(real(ifft(u)).*ifft(conj(k).*v));
+
+end

src/+otp/+kuramotosivashinsky/jvp.m

@@ -0,0 +1,5 @@
+function jv = jvp(~, u, k, k2, k4, v)
+
+jv = -k2.*v - k4.*v - k.*fft(real(ifft(u)).*ifft(v));
+
+end