问题描述
我正在尝试找到一种使此代码更快的方法。
I am trying to find a way to make this code faster.
Nagumo1是用于计算时间t处两个导数的值的函数。
Nagumo1 is the function that calculates the value of the two derivatives at time 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）以毫秒为单位。
It is a system of differential equations that represents a largely simplified model of neuron. V represents a difference of electric potential, n represents the number of K+/Na+ canals and Iapp is the electric current applied to the neuron. The time variable (t) is measured in msec.
我想使用具有可变步长的Euler显式方法来数值求解微分方程组和图形
I want to use the Euler explicit method, with a variable step size, to numerically resolve the differential equation system and graphe the solution.
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解决方案相比（我认为是准确的）相差很大。
I haven't included any comments, because I think it is clear. I just want to maintain the adaptable step h and make the code faster. Ideally I would like to find a way to initialize z and temps(time), but when I try to do that then I have a problem plotting my solution. Note that when erreura(absolute error) and erreurr(relative error) are greater than 10^-6 my solution varies a lot in comparison to ode45 solution (which i consider to be accurate).
有什么想法吗？
PS如果要测试使用的值在-2，V之间为2、2，n为0.1、1，Iapp为1（定义了一个函数句柄@（t））之间变化，则为-2。
P.S. if you want to test use values varying between -2, 2 for V, 0,1, 1 for n, 0.1, 1 for Iapp (defined a function handle @(t)).
推荐答案
在尝试加快解释后的代码的速度之前，您应该注意获得一个正确的解决方案。在 to + i * h 的时间计算中仍然存在一些不对的地方，仅对固定步长有效。我将从基本原理开始解释自适应方法。
Before trying to speed up the interpreted code, you should care to get a correct solution at all. That there is still something amiss is visible in the time computations to+i*h that are only valid for a fixed step size. I'll explain the adaptive method from first principles.
使用的近似值是时间 t 以步长 h 计算得出的一阶精确解与
uses the approximation that the numerical solution at time t computed with step size h relates to the exact solution in first order as
y(h;t)=y_exact(t) + C*t*h + O(t*h²)
认为前进一到两步的一半有错误
gives that the advancement in one and two steps of half size has the errors
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 即可得到与全局误差直接比较的值。
is an estimator of the local error at step size h/2. We know that the local errors in first order add to the global error (in a better approximation there is some compounding with the Lipschitz constant as "annual" interest rate). Thus in the reverse direction we desire to get that the local error is a h sized part of the global error. Divide all local error quantities by h to get values that directly compare to the global error.
现在尝试保持该局部误差估计 local_err = norm（y（h; h）-y（h / 2; h）） / h =某个走廊内的标准（C）* h / 2 [tol / 100，tol] 其中 tol代表所需的全局变量错误。因此，根据当前数据计算出的理想步长为
now tries to keep that local error estimate local_err = norm(y(h;h)-y(h/2;h))/h = norm(C)*h/2 inside some corridor [tol/100, tol] where ´tol´ stands for the desired global error. The ideal step size from the current data is thus computed as
tol = norm(C)*h_ideal/2 = local_err*h_ideal/h
<==>
h_ideal = tol / local_err * h
在算法中，这些积分步骤和误差估计，然后接受这些步骤，如果在公差范围内，则继续进行计算，然后根据上述公式调整步长，以进行循环的下一次迭代。除了使用计算出的理想步长，还可以通过恒定因子沿理想步长的方向修改步长。通常，这只会增加被拒绝的步数，从而仍然达到理想的步长。
In the algorithm one would compute these integration steps and error estimates and then accept the steps and advance the computation if inside the tolerance bounds, then adapt the step size by above formula go to the next iteration of the loop. Instead of using the computed ideal step size one could also modify the step size by constant factors in the direction of the ideal step size. In general this will only increase the number of rejected steps to still reach the ideal step size.
为避免尝试的和使用的步长出现振荡和过于突然的变化，请引入某种移动平均值，抑制方向 1 的变化因子，例如
To avoid oscillations and too abrupt changes in the tried and used step sizes, introduce some kind of moving average, dampen the change factor in direction 1 like in
a = tol / (1e-12+local_err);
h = 0.25*(3+a^0.8)*h ;
在代码中，这看起来像是
In code this could look like
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
此方法的实际应用给出了以下图表。
The practical application of this method gives the following plot.
在顶部显示了求解曲线。人们看到在弯曲的或快速变化的部分上有较高的密度，而在溶液曲线更平直的地方则看到较低的密度。在下部显示针对最小公差的解决方案的误差。差异是根据解决方案的容忍度来缩放的，因此所有这些都共享相同的缩放比例。可以看到，输出紧密跟踪输入所要求的公差。
In the top the solution curves are depicted. One sees a higher density at the curved or rapidly changing parts and a lower density where the solution curve is more straight. In the lower part the error against the solution of lowest tolerance is displayed. The difference is scaled by the tolerance of the solution, so that all share the same scale. As one can see, the output traces the tolerance demanded at input closely.
这篇关于具有自适应步骤的Matlab Euler Explicit ode求解器，有没有办法使代码更快？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！