具有自适应步长的 Matlab Euler Explicit ode 求解器，有没有办法使代码更快？
Matlab Euler Explicit ode solver with adaptable step, is there a way to make code faster?
我正试图找到一种方法来加快这段代码的速度。
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+ 通道的数量，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 解决方案（我认为是准确的）相比变化很大。
有什么想法吗？
P.S。如果你想测试使用值在 -2、V 为 2、n 为 0,1、1、Iapp 为 0.1、1 之间变化（定义函数句柄 @(t)）。
看看那些行：
z = [ z ; y ];
temps = [temps; temps(i)+h];
这些真的很慢，而且我知道在使用可变步长时你不能 pre-allocate。假设您使用的数据量很大，我建议您使用传统文件。 z 的替代品是：
fp = fopen('z_file.dat','wb');
...
fwrite(fp,y,'double');
...
fclose(fp);
fp = fopen('z_file.dat','r');
z = fread(fp,length_of_z,'double');
fclose(fp);
您需要知道存储数据量的地方（length_of_z 我猜是您的 "i"）。此外，如果数据量相当大，这只会加快速度。
在尝试加速解释代码之前，您应该注意获得正确的解决方案。在仅对固定步长有效的时间计算 to+i*h 中仍然存在一些问题。我将从第一原则解释自适应方法。
理查森外推法的误差估计
使用步长 h 在时间 t 计算的数值解与一阶精确解相关的近似值作为
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 常数的复合，如 "annual" 利率）。因此，在相反的方向上，我们希望得到局部误差是全局误差的 h 大小的一部分。将所有局部误差量除以 h 以获得与全局误差直接比较的值。
自适应步长控制器
现在尝试将局部误差估计 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
该方法的实际应用给出如下图
在顶部描绘了解曲线。在弯曲或快速变化的部分可以看到更高的密度，而在解曲线更直的地方可以看到更低的密度。在下部显示了针对最低公差解决方案的错误。差异由解决方案的公差缩放，以便所有共享相同的比例。可以看出，输出密切跟踪输入要求的公差。
我正试图找到一种方法来加快这段代码的速度。
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+ 通道的数量，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 解决方案（我认为是准确的）相比变化很大。
有什么想法吗？
P.S。如果你想测试使用值在 -2、V 为 2、n 为 0,1、1、Iapp 为 0.1、1 之间变化（定义函数句柄 @(t)）。
看看那些行：
z = [ z ; y ];
temps = [temps; temps(i)+h];
这些真的很慢，而且我知道在使用可变步长时你不能 pre-allocate。假设您使用的数据量很大，我建议您使用传统文件。 z 的替代品是：
fp = fopen('z_file.dat','wb');
...
fwrite(fp,y,'double');
...
fclose(fp);
fp = fopen('z_file.dat','r');
z = fread(fp,length_of_z,'double');
fclose(fp);
您需要知道存储数据量的地方（length_of_z 我猜是您的 "i"）。此外，如果数据量相当大，这只会加快速度。
在尝试加速解释代码之前，您应该注意获得正确的解决方案。在仅对固定步长有效的时间计算 to+i*h 中仍然存在一些问题。我将从第一原则解释自适应方法。
理查森外推法的误差估计
使用步长 h 在时间 t 计算的数值解与一阶精确解相关的近似值作为
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 常数的复合，如 "annual" 利率）。因此，在相反的方向上，我们希望得到局部误差是全局误差的 h 大小的一部分。将所有局部误差量除以 h 以获得与全局误差直接比较的值。
自适应步长控制器
现在尝试将局部误差估计 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
该方法的实际应用给出如下图
在顶部描绘了解曲线。在弯曲或快速变化的部分可以看到更高的密度，而在解曲线更直的地方可以看到更低的密度。在下部显示了针对最低公差解决方案的错误。差异由解决方案的公差缩放，以便所有共享相同的比例。可以看出，输出密切跟踪输入要求的公差。