Improvements for Jacobian and JAVP

Author: Steven-Roberts

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

@@ -7,8 +7,8 @@
 % 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." 
+% 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.
 %
 
@@ -20,26 +20,22 @@
             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;
-            
             % exclude the left boundary point as it is identical to the
             % right boundary point
-            x = linspace(h, L, N).';
+            x = linspace(L / N, L, N).';
 
-            u0 = cospi(2*x/L).*(1+sinpi(2*x/L));
+            u0 = cospi(2 * x / L) .* (1 + sinpi(2 * x / L));
 
             u0hat = fft(u0);
-
+            
             tspan = [0, 150];
 
             obj = [email protected](tspan, u0hat, params);

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

@@ -5,8 +5,8 @@
 
             p = inputParser;
             addParameter(p, 'Size', 1024);
-            addParameter(p, 'L', 80*pi);
-
+            addParameter(p, 'L', 80 * pi);
+            
             parse(p, varargin{:});
 
             s = p.Results;
@@ -16,21 +16,19 @@
 
             params.L = L;
 
-            h = L/N;
-            
             % exclude the left boundary point as it is identical to the
             % right boundary point
-            x = linspace(h, L, N).';
+            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, max(eta2, max(eta3, eta4))));
+            u0 = 16*max(0, max([eta1, eta2, eta3, eta4], [], 2));
 
             u0hat = fft(u0);
-
+            
             tspan = [0, 100];
 
             obj = [email protected](tspan, u0hat, params);

src/+otp/+kuramotosivashinsky/KuramotoSivashinskyProblem.m

@@ -9,35 +9,27 @@
     methods
         function soly = convert2grid(~, soly)
 
-            soly = real(ifft(soly));
+            soly = abs(ifft(soly));
 
         end
 
     end
 
     methods (Access = protected)
         function onSettingsChanged(obj)
-            
-            L = obj.Parameters.L;
-            
             N = obj.NumVars;
 
-            div = L/(2*pi);
+            k = 2 * pi * [0:(N/2 - 1), 0, (-N/2 + 1):-1].' / obj.Parameters.L;
+            ik = 1i * k;
+            k24 = k.^2 - k.^4;
 
-            % note that k already has "i" in it
-            k = (1i*[0:(N/2 - 1), 0, (-N/2 + 1):-1].'/div);
-            k2 = k.^2;
-            k4 = k.^4;
-            k24 = k2 + k4;
-
-            obj.Rhs = otp.Rhs(@(t, u) otp.kuramotosivashinsky.f(t, u, k, k24), ...
+            obj.Rhs = otp.Rhs(@(t, u) otp.kuramotosivashinsky.f(t, u, ik, k24), ...
                 otp.Rhs.FieldNames.Jacobian, ...
-                @(t, u) otp.kuramotosivashinsky.jac(t,u, k, k24), ...
+                @(t, u) otp.kuramotosivashinsky.jac(t,u, ik, k24), ...
                 otp.Rhs.FieldNames.JacobianVectorProduct, ...
-                @(t, u, v) otp.kuramotosivashinsky.jvp(t, u, v, k, k24), ...
+                @(t, u, v) otp.kuramotosivashinsky.jvp(t, u, v, ik, k24), ...
                 otp.Rhs.FieldNames.JacobianAdjointVectorProduct, ...
-                @(t, u, v) otp.kuramotosivashinsky.javp(t, u, v, k, k24));
-            
+                @(t, u, v) otp.kuramotosivashinsky.javp(t, u, v, ik, k24));
         end
 
         function validateNewState(obj, newTimeSpan, newY0, newParameters)

src/+otp/+kuramotosivashinsky/f.m

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

src/+otp/+kuramotosivashinsky/jac.m

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

src/+otp/+kuramotosivashinsky/javp.m

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

src/+otp/+kuramotosivashinsky/jvp.m

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