Implement 2N low-storage methods.

Author: ketch

RK-coeff-opt/nonlinear_constraints.m

@@ -6,7 +6,7 @@
 % The input arguments are:
 %     * :math:`x`: vector of the decision variables.  See unpack_rk.m for details about
 %       the order in which they are stored.
-%     * *class*: class of method to search ('erk' = explicit RK; 'irk' = implicit RK; 'dirk' = diagonally implicit RK; 'sdirk' = singly diagonally implicit RK; '2S', '3S', '2S*', '3S*' = low-storage formulations).
+%     * *class*: class of method to search ('erk' = explicit RK; 'irk' = implicit RK; 'dirk' = diagonally implicit RK; 'sdirk' = singly diagonally implicit RK; '2S', '3S', '2S*', '3S*', '2N' = low-storage formulations).
 %     * :math:`s`:number of stages.
 %     * :math:`p`: order of the RK scheme.
 %     * *objective*: objective function ('ssp' = maximize SSP coefficient; 'acc' = minimize leading truncation error coefficient).

RK-coeff-opt/rk_obj.m

@@ -5,7 +5,7 @@
 %
 % Inputs:
 %     * :math:`x`: vector of the unknowns.
-%     * class: class of method to search ('erk' = explicit RK; 'irk' = implicit RK; 'dirk' = diagonally implicit RK; 'sdirk' = singly diagonally implicit RK; '2S', '3S', '2S*', '3S*' = low-storage formulations).
+%     * *class*: class of method to search ('erk' = explicit RK; 'irk' = implicit RK; 'dirk' = diagonally implicit RK; 'sdirk' = singly diagonally implicit RK; '2S', '3S', '2S*', '3S*', '2N' = low-storage formulations).
 %     * :math:`s`:number of stages.
 %     * :math:`p`: order of the RK scheme.
 %     * objective: objective function ('ssp' = maximize SSP coefficient; 'acc' = minimize leading truncation error coefficient).

RK-coeff-opt/rk_opt.m

@@ -150,7 +150,7 @@ pctRunOnAll warning('on', 'all')
 
 
 % Check order of the scheme
-if strcmp(class(1:2),'2S') || strcmp(class(1:2),'3S')
+if strcmp(class(1:2),'2S') || strcmp(class(1:2),'3S' || strcmp(class), '2N')
     [rk.A,rk.Ahat,rk.b,rk.bhat,rk.c,rk.chat,rk.alpha,rk.beta,rk.gamma1,rk.gamma2,rk.gamma3,rk.delta] = unpack_lsrk(X,class);
     order = check_RK_order(rk.A,rk.b,rk.c,'nonlinear');
 elseif k==1

RK-coeff-opt/set_n.m

@@ -34,6 +34,8 @@
       n = 4*s - 3;
     case '3SstarembFSAL'      %Low-storage 3S* embedded pairs of Ketcheson
       n = 4*s - 2;
+    case '2N'             % Williamson
+      n = 2*s -1;
 
     %=====================
     % Multistep RK classes

RK-coeff-opt/unpack_lsrk.m

@@ -50,6 +50,26 @@
     %     gamma_{53} ... gamma_{s+1,3} delta_s delta_{s+1} delta_{s+2} ]
     delta =[1 X(s:2*s-3) X(end-3:end-1)];
     betahat = X(end);
+  case '2N'
+    % n = 2s - 1 free parameters
+    s = (length(X)+1)/2; % # of stages
+    % for 2N methods:
+    % X = [A_2, A_3, ..., A_s, B_1, ..., B_s]
+    A_2N = [0 X(1:s-1)];
+    B_2N = X(s:end);
+    K = zeros(s+1,s);
+    for i=1:s-1
+        K(i,i+1) = B_2N(i);
+        for j = 2:s-i+1
+            for k=0:j-1
+                K(i+j,i) = K(i+j,i) + B_2N(i+k)*prod(A_2N(i+1:i+k));
+            end
+        end
+    end
+    A = K(1:s,:);
+    b = K(s+1,:); b=b';
+    c=sum(A,2);
+
 end
 
 if strcmp(class(1:2),'2S')

RK-coeff-opt/unpack_rk.m

@@ -88,6 +88,26 @@
         %     gamma_{53} ... gamma_{s+1,3} delta_s delta_{s+1} delta_{s+2} ]
         delta =[1 X(s:2*s-3) X(end-3:end-1)];
         betahat = X(end);
+    case '2N'
+        % n = 2s - 1 free parameters
+        s = (length(X)+1)/2; % # of stages
+        % for 2N methods:
+        % X = [A_2, A_3, ..., A_s, B_1, ..., B_s]
+        A_2N = [0 X(1:s-1)];
+        B_2N = X(s:end);
+        K = zeros(s+1,s);
+        for i=1:s-1
+            K(i,i+1) = B_2N(i);
+            for j = 2:s-i+1
+                for k=0:j-1
+                    K(i+j,i) = K(i+j,i) + B_2N(i+k)*prod(A_2N(i+1:i+k));
+                end
+            end
+        end
+        A = K(1:s,:);
+        b = K(s+1,:); b=b';
+        c=sum(A,2);
+
 end