Confirmation of FSAL property for IMEX methods by Kennedy and Carpenter

Question

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    

Answer 1

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.