This question is a continuation of Fourth order IMEX Runge-Kutta method and Implementation details for high order IMEX methods by Kennedy and Carpenter.
I need confirmation that ARK3(2)4L[2]SA by Kennedy and Carpenter has the "First Same As Last"-property. I have verified it numerically for ARK4(3)6L[2]SA and ARK5(4)8L[2]SA, but do not get it for ARK3(2)4L[2]SA. I have attached my MatLab-code, and if you run in you will see that differenceStep in function ARK3 is not zero.
function [] = testARK()
a = -1;
b = 1;
c = 2;
f_s = @(t,y) a*c*y;
f_ns = @(t,y) 2*b*t;
f = @(t,y) f_s(t,y) + f_ns(t,y);
solution = @(t) c*exp(a*c*t) + (2*b*(-1 + exp(a*c*t) - a*c*t))/(a^2*c^2);
% Initial data
t0 = 0;
T = 1; % Final time
y0 = solution(t0);
dt_v = T*2.^(-1:-1:-9);
err_V = zeros(length(dt_v),1);
err_est_V = zeros(length(dt_v),1);
for i = 1:length(dt_v)
yn = y0;
t = 0;
dt = dt_v(i);
ks_4 = f_s(t0,y0);
while t<T
[yn,err_est, t,dt,ks_4] = ARK3(yn,t,dt,f_s, f_ns,ks_4);
end
err_V(i) = abs(yn-solution(t));
err_est_V(i) = err_est;
end
figure(1)
loglog(dt_v,err_V,'-')
hold on
loglog(dt_v,err_est_V,':')
loglog(dt_v,dt_v.^(3),'o-')
hold off
leg = legend('err','err estimate','$\mathcal{O}(\Delta t^{3})$');
set(leg,'interpreter','latex');
title('ARK3')
Table = table(err_V(1:end-1)./err_V(2:end));
Table.Properties.VariableNames = {'RateOfDecay'};
disp(Table)
end
function [yn,err_est,t,dt,ks_4] = ARK3(yn,t,dt,f_s, f_ns,ks_4)
%% ARK3(2)4L[2]SA-ERK
A_E = [0,0,0,0;...
1767732205903/2027836641118, 0 ,0, 0;...
5535828885825/10492691773637,788022342437./10882634858940, 0, 0;...
6485989280629/16251701735622,-4246266847089/9704473918619,...
10755448449292/10357097424841,0];
b = [1471266399579/7840856788654,...
-4482444167858/7529755066697, ...
11266239266428/11593286722821,...
1767732205903/4055673282236];
b_hat = [2756255671327/12835298489170,...
-10771552573575/22201958757719,...
9247589265047/10645013368117,...
2193209047091/5459859503100];
c = [0;1767732205903/2027836641118;3/5;1];
%% ARK4(3)6L[2]SA-ESDIRK
A_I = [0,0,0,0;...
1767732205903/4055673282236,1767732205903/4055673282236,0,0;...
2746238789719/10658868560708,-640167445237/6845629431997,...
1767732205903/4055673282236,0;...
b];
%%
s = 4; % Number of stages
k_v = zeros(s,2); % First column: stiff. Second column: non-stiff
k_v(1,1) = f_s(t,yn);
differenceStep = abs(f_s(t,yn)-ks_4)
k_v(1,2) = f_ns(t,yn);
for i = 2:s
k_v(i,1) = f_s(t + dt * c(i), yn +...
dt * A_E(i,1:(i-1))*k_v(1:(i-1),2) +...
dt * A_I(i,1:(i-1))*k_v(1:(i-1),1))...
/(1-f_s(t + dt * c(i), dt*A_I(i,i)));
k_v(i,2) = f_ns(t + dt * c(i), yn +...
dt * A_E(i,1:(i-1))*k_v(1:(i-1),2) +...
dt * A_I(i,1:i)*k_v(1:i,1));
end
ks_4 = k_v(s,1);
yn_hat = yn + dt*b_hat*sum(k_v,2);
yn = yn + dt*b*sum(k_v,2);
err_est = abs(yn-yn_hat);
t = t + dt;
end