# Confirmation of FSAL property for IMEX methods by Kennedy and Carpenter

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


## 1 Answer

I don't think I was clear so let me clarify my previous post. The FSAL property is on the implicit tableau of the ESDIRK method. When used as an ESDIRK method it is FSAL. When the IMEX portion is in there, it's not ESDIRK because the updated dependent variable adds in the explicit portion.

• That makes sense. So FSAL actually can't be used for IMEX methods?. However, it works perfectly for ARk4 and ARK5, why is that? – Raibyo Mar 14 '18 at 20:27
• I honestly don't know. – Chris Rackauckas Mar 15 '18 at 0:37