0
$\begingroup$

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    
$\endgroup$
0
$\begingroup$

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.