我正在尝试找到一种方法来使此代码更快。
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 解决方案（我认为是准确的）相比变化很大。
有任何想法吗？
附：如果你想测试使用-2、2（V）、0,1、1（n）、0.1、1（Iapp）之间的变化值（定义了一个函数句柄@(t)）。
在尝试加速解释代码之前，您应该注意获得正确的解决方案。从时间计算中可以看出仍然存在一些问题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³)
and thus
y(h;h)-y(h/2;h) = C*h²/2 + O(h³)
是步长下局部误差的估计h/2。我们知道，一阶局部误差会增加全局误差（在更好的近似中，存在与利普希茨常数作为“年”利率的复合）。因此，在相反的方向上，我们希望得到局部误差是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
该方法的实际应用如下图所示。
顶部描绘了解决方案曲线。人们会在弯曲或快速变化的部分看到较高的密度，而在解曲线更直的地方看到较低的密度。在下部显示相对于最低公差解的误差。差异按解决方案的容差进行缩放，以便所有的都具有相同的比例。正如我们所看到的，输出密切跟踪输入所需的容差。