<p>最后一点:</p>
<ul>
<li><p>在当前表单中,它甚至不是
端庄的颂歌。它过早地停止了一个步骤,即最后一个“常规”步骤
应该朝向<code>noiter*dt</code>,并且不考虑时间
余数<code>T-noiter*dt</code>。在</p>
<p>注意,<code>range(N)</code>生成的是<code>0,1,…,N-1</code>。同样地,
<code>res=zeros(N)</code>生成一个包含<code>N</code>项的数组,从<code>res[0]</code>到
<code>res[N-1]</code>。</p></li>
<li><p>转换不应依赖于离散化,即步骤
长度。为了达到这个目的,一个更精确的穿越时间
应通过插值(线性或
反向二次方),然后用
新的初始条件。要保留所需的栅格,请使用short
第一步。</p></li>
</ul>
<hr/>
<pre><code>def solve_my_ODE(res, fdot, x0, T, dt):
noiter = int(T / dt)
dt = T/noiter #adapt the timestep
res = zeros(noiter+1)
res[0] = x0
for i in range(noiter):
res[i + 1] = res[i] + dt * fdot(res[i], i * dt)
if res[i + 1] >= 2:
h = (2-res[i])/(res[i+1]-res[i]) # precautions against zero division ?
res[i + 1] = 0 + (1-h)*dt * fdot(0, (i+h)*dt)
return res
</code></pre>
<hr/>
<ul>
<li><p>看来,最终精度要优于<code>1e-4</code>。
这里使用<code>dt=1e-5</code>计算使用<code>100 000</code>步骤
以及同样多的功能评估。在</p>
<p>使用经典的Runge-Kutta方法和<code>h=0.05</code>will
导致错误略大于<code>1e-5</code>(<code>dt**4=6.25e-6</code>),
i、 与欧拉法的误差大小相当。
但是,现在这只需要<code>T/dt=20</code>个步骤,总共<code>80</code>
功能评估。注意,切换时间也需要
以<code>O(dt**4)</code>顺序精确,以防污染
全局错误顺序。在</p>
<p>因此,如果以速度为目标,研究
高阶方法。</p></li>
</ul>