Kennedy 和 Carpenter 确认 IMEX 方法的 FSAL 属性

计算科学 时间积分 龙格库塔 自适应时间步长
2021-12-09 02:29:56

这个问题是 Kennedy 和 Carpenter 的四阶 IMEX Runge-Kutta 方法高阶 IMEX 方法的实施细节的延续

我需要确认 Kennedy 和 Carpenter 的 ARK3(2)4L[2]SA 具有“First As Last”属性。我已经对 ARK4(3)6L[2]SA 和 ARK5(4)8L[2]SA 进行了数值验证,但对于 ARK3(2)4L[2]SA,没有得到它。我附上了我的 MatLab 代码,如果你在其中运行,你会看到函数 ARK3 中的 differenceStep 不为零。

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个回答

我不认为我很清楚,所以让我澄清一下我之前的帖子。FSAL 属性位于 ESDIRK 方法的隐式画面上。当用作 ESDIRK 方法时,它是 FSAL。当 IMEX 部分在那里时,它不是 ESDIRK,因为更新的因变量添加到显式部分中。