Completely changed the problem to be spectral

Author: AndreyAPopov

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

@@ -4,26 +4,30 @@
         function obj = Canonical(varargin)
 
             p = inputParser;
-            addParameter(p, 'Size', 64, @isscalar);
-            addParameter(p, 'L', 25, @isscalar);
+            addParameter(p, 'Size', 128);
+            addParameter(p, 'L', 32*pi);
 
             parse(p, varargin{:});
 
             s = p.Results;
 
-            n = s.Size;
+            N = s.Size;
+            L = s.L;
 
+            params.L = L;
 
-            params.n = n;
-            params.l = s.L;
+            h=L/N;
 
-            x = linspace(0, 10*pi, n + 1);
+            x=h*(1:N).';
 
-            u0 = 4*cos(x(1:end-1)).';
+            u0 = cos(x/16).*(1+sin(x/16));
+            
+            u0hat = fft(u0);
 
             tspan = [0, 100];
 
-            obj = obj@otp.kuramotosivashinsky.KuramotoSivashinskyProblem(tspan, u0, params);
+            obj = obj@otp.kuramotosivashinsky.KuramotoSivashinskyProblem(tspan, u0hat, params);
+            
         end
     end
 end

src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m

@@ -2,31 +2,57 @@
 
     methods
         function obj = KuramotoSivashinskyProblem(timeSpan, y0, parameters)
-            obj@otp.Problem('Kuramoto Sivashinsky Problem', [], timeSpan, y0, parameters);
+            obj@otp.Problem('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)
-            n = obj.Parameters.n;
-            l = obj.Parameters.l;
 
-            domain = [-l, l];
-            D = otp.utils.pde.D(n, domain, 'C');
-            L = otp.utils.pde.laplacian(n, domain, 1, 'C');
+            L = obj.Parameters.L;
 
-            obj.Rhs = otp.Rhs(@(t, u) otp.kuramotosivashinsky.f(t, u, D, L), ...
-                otp.Rhs.FieldNames.Jacobian, @(t, u) otp.kuramotosivashinsky.jac(t, u, D, L));
+            N = numel(obj.Y0);
+            
+            div = L/(2*pi);
+            
+            k = (1i*[0:(N/2 - 1), 0, (-N/2 + 1):-1].'/div);
+            k2 = k.^2;
+            k4 = k.^4;
+            
+            Dfft = dftmtx(N);
+            Difft = conj(Dfft)/N;
+
+            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, Dfft, Difft), ...
+                otp.Rhs.FieldNames.JacobianVectorProduct, ...
+                @(t, u, v) otp.kuramotosivashinsky.jvp(t,u, k, k2, k4, v));
 
         end
 
         function validateNewState(obj, newTimeSpan, newY0, newParameters)
+            
             validateNewState@otp.Problem(obj, ...
                 newTimeSpan, newY0, newParameters)
 
+            assert(mod(numel(newY0), 2) == 0);
+            
             otp.utils.StructParser(newParameters) ...
-                .checkField('n', 'scalar', 'integer', 'finite', 'positive') ...
-                .checkField('l', 'scalar', 'integer', 'finite', 'positive');
+                .checkField('L', 'scalar', 'finite', 'positive');
+            
         end
     end
 end

src/+otp/+kuramotosivashinsky/f.m

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

src/+otp/+kuramotosivashinsky/jac.m

@@ -1,5 +1,5 @@
-function j = jac(~, u, D, L)
+function j = jac(~, u, k, k2, k4, Dfft, Difft)
 
-j = -L*L - L - spdiags(u, 0, numel(u), numel(u))*D - spdiags(D*u, 0, numel(u), numel(u));
+j = -diag(k2) - diag(k4) - diag(k)*Dfft*diag(ifft(u))*Difft;
 
 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)).*real(ifft(v)));
+
+end