differential-equations - 具有自适应步长的 Matlab Euler Explicit ode 求解器，有没有办法使代码更快？

Question

我正试图找到一种方法来使这段代码更快。

Nagumo1 是计算两个导数在时间 t 的值的函数。

function x = nagumo(t, y, f)

Iapp = f(t);
e = 0.1;
F = 2/(1+exp(-5*y(1)));
n0 = 0;

x = zeros(2, 1);

z(1) = y(1) - (y(1).^3)/3 - y(2).^2 + Iapp;  %z(1) = dV/dt

z(2) = e.*(F + n0 - y(2));                   %z(2) = dn/dt

x = [z(1);z(2)];
end

它是一个微分方程组，代表了一个大大简化的神经元模型。 V 表示电位差，n 表示 K+/Na+ channel 的数量，Iapp 是施加到神经元的电流。时间变量 (t) 以毫秒为单位。

我想使用可变步长的欧拉显式方法对微分方程组进行数值求解并绘制解。

function x =  EulerExplicit1(V0, n0, tspan, Iapp) 
 format long;

 erreura = 10^-3;
 erreurr = 10^-6;

 h = 0.1;                             

 to =tspan(1, 1) + h;                 
 temps = tspan(1, 1);
 tf = tspan(1, 2);

 y = zeros(1,2);
 yt1 = zeros(1, 2);
 yt2 = zeros(1, 2);
 y = [V0, n0];  

 z = y;

 i = 1;

 s = zeros(1, 2);
 st1 = zeros(1, 2);

 while temps<tf

     s = nagumo1(to+i*h, y, Iapp);
     y = y + h*s;
     yt1 = y + (h/2)*s;
     st1 = nagumo1(to+(i*h+h/2), yt1, Iapp);
     yt2 = yt1 + (h/2)*st1;

     if abs(yt2-y)>(erreura+erreurr*abs(y))
        test = 0;
     elseif h<0.4
         h = h*2;
         test = 0;
     end

     while test == 0

         if abs(yt2-y)>(erreura+erreurr*abs(y)) & h>0.01
             h = h/2;
             s = nagumo1(to+i*h, y, Iapp);
             y = y + h*s;
             yt1 = y + (h/2)*s;
             st1 = nagumo1(to+i*h+h/2, yt1, Iapp);
             yt2 = yt1 + (h/2)*st1;
         else
             test = 1;
         end
     end
     z = [ z ; y ];
     temps = [temps; temps(i)+h];
     i = i+1;

 end

 x = zeros(size(z));
 x = z;

 disp('Nombre d iterations:');
 disp(i);
 plot(temps, x(:, 1:end), 'y');
 grid;

end

我没有包含任何评论，因为我认为它很清楚。我只想保持可适应的步骤 h 并使代码更快。理想情况下，我想找到一种方法来初始化 z 和 temps(time)，但是当我尝试这样做时，我在绘制解决方案时遇到了问题。请注意，当 erreura(绝对误差)和 erreurr(相对误差)大于 10^-6 时，我的解决方案与 ode45 解决方案(我认为是准确的)相比变化很大。

有什么想法吗？

附言如果你想测试使用值在 -2、V 为 2、n 为 0,1、1、Iapp 为 0.1、1 之间变化(定义函数句柄 @(t))。

Answer 1

在尝试加速解释代码之前，您应该注意获得正确的解决方案。在仅对固定步长有效的时间计算 to+i*h 中仍然存在一些错误是可见的。我将从第一原则解释自适应方法。

理查森外推法的误差估计

使用时间 t 以步长 h 计算的数值解与一阶精确解相关的近似值作为

y(h;t)=y_exact(t) + C*t*h + O(t*h²)

给出半尺寸的一步和两步的进步有错误

y(h;h) = y_exact(h) + C*h² + O(h³)
y(h/2;h) = y_exact(h)+C*h²/2 + O(h³)

因此

y(h;h)-y(h/2;h) = C*h²/2 + O(h³)

是步长 h/2 的局部误差估计量。我们知道一阶局部误差会添加到全局误差中(在更好的近似中，有一些与 Lipschitz 常数的复合作为“年”利率)。因此，在相反的方向上，我们希望得到局部错误是全局错误的 h 大小的一部分。将所有局部误差量除以 h 以获得与全局误差直接比较的值。

自适应步长 Controller

现在尝试保留局部误差估计 local_err = norm(y(h;h)-y(h/2;h))/h = norm(C)*h/2一些走廊 [tol/100, tol] 其中“tol”代表所需的全局误差。因此，当前数据的理想步长计算为

 tol = norm(C)*h_ideal/2 = local_err*h_ideal/h

 <==>

 h_ideal = tol / local_err * h 

在算法中，人们会计算这些积分步骤和误差估计，然后接受步骤并在公差范围内推进计算，然后通过上述公式调整步长以进入循环的下一次迭代。除了使用计算出的理想步长之外，还可以在理想步长的方向上通过常数因子修改步长。一般来说，这只会增加被拒绝的步数，以达到理想的步长。

为避免尝试和使用的步长出现振荡和过于突然的变化，引入某种移动平均线，抑制方向1的变化因子，如

 a = tol / (1e-12+local_err);
 h = 0.25*(3+a^0.8)*h ;

在代码中这看起来像

while t < t_final
    if t+1.1*h > t_final
        h = t_final - t
        force_advance = True
    end 
    s1  = f(t,y)
    s05 = f(t+0.5*h, y+0.5*h*s1)
    s2  = 0.5*(s1+s05)

    localerr = norm(s2-s1)
    tol = abstol+norm(y)*reltol
    if force_advance | (0.01*tol < localerr & localerr < tol)
        y = y + h*s2
        t = t + h
        sol_y(end+1)=y
        sol_t(end+1)=t
        force_advance = False
    end
    a = tol / (1e-19+localerr) )
    h = 0.25*(3+a^0.8)*h ;
    if h < h_min
        h = h_min
        force_advance = True
    end
    if h > h_max
        h = h_max
        force_advance = True
    end
end

该方法的实际应用给出了如下图。

在顶部描绘了解曲线。在弯曲或快速变化的部分可以看到更高的密度，而在解曲线更直的地方可以看到更低的密度。在下部显示了针对最低公差解决方案的错误。差异由解决方案的公差缩放，以便所有共享相同的比例。可以看出，输出密切跟踪输入要求的公差。