<h3>物理</h3>
<p>牛顿定律给你一个二阶常微分方程<code>u''=F(u)</code>,其中<code>u=[x,y]</code>。使用<code>v=[x',y']</code>得到一阶系统</p>
<pre><code>u' = v
v' = F(u)
</code></pre>
<p>它是四维的,必须用四维状态来求解。唯一可用的简化方法是使用开普勒定律,它允许将系统缩减到一个标量级的角度。但这不是这里的任务。在</p>
<h3>欧拉法</h3>
<p>您正确地实现了Euler方法来计算代码的最后一个循环中的值。它看起来不是物理的可能是因为Euler方法不断地增加轨道,因为它沿着切线移动到凸轨迹的外部。在您的实现中,可以看到<code>G=100</code>的这种向外螺旋。在</p>
<p><a href="https://i.stack.imgur.com/DZ3Fp.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/DZ3Fp.png" alt="enter image description here"/></a></p>
<p>通过选择较小的步长,例如<code>dt=0.001</code>,可以有效地减少这种情况。在</p>
<p><a href="https://i.stack.imgur.com/o8KEz.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/o8KEz.png" alt="enter image description here"/></a></p>
<p>你应该选择积分时间作为一个完整轨道的一个很好的部分,以得到一个体面的结果,与上述参数,你得到大约2个循环,这是好的。在</p>
<h3>RK4实施</h3>
<p>你犯了几个错误。不知怎么的你失去了速度,位置更新应该基于速度。在</p>
<p>那么您应该在<code>fx(x + .5*kx1, y + .5*kx1, t + .5*dt)</code>处停下来重新考虑您的方法,因为这与任何命名约定都不一致。一致、正确的变体是</p>
^{pr2}$
<p>这表明您不能解耦耦合系统的集成,因为您需要在<code>x</code>更新的同时进行<code>y</code>更新。此外,函数值是加速度,从而更新速度。位置更新使用当前状态的速度。因此,步骤应该从</p>
<pre><code> kx1 = dt * fx(x,y,t) # vx update
mx1 = dt * vx # x update
ky1 = dt * fy(x,y,t) # vy update
my1 = dt * vy # y update
kx2 = dt * fx(x + 0.5*mx1, y + 0.5*my1, t + 0.5*dt)
mx2 = dt * (vx + 0.5*kx1/2)
ky2 = dt * fy(x + 0.5*mx1, y + 0.5*my1, t + 0.5*dt)
my2 = dt * (vy + 0.5*ky1/2)
</code></pre>
<p>等等</p>
<p>然而,正如您所见,这已经开始变得难以处理。将状态组合成一个向量,并对系统方程使用向量值函数</p>
<pre><code>M, G = 20, 100
def orbitsys(u):
x,y,vx,vy = u
r = np.hypot(x,y)
f = G*M/r**3
return np.array([vx, vy, -f*x, -f*y]);
</code></pre>
<p>然后可以使用Euler或Runge-Kutta步骤的烹饪书实现</p>
<pre><code>def Eulerstep(f,u,dt): return u+dt*f(u)
def RK4step(f,u,dt):
k1 = dt*f(u)
k2 = dt*f(u+0.5*k1)
k3 = dt*f(u+0.5*k2)
k4 = dt*f(u+k3)
return u + (k1+2*k2+2*k3+k4)/6
</code></pre>
<p>并将它们组合成一个集成回路</p>
<pre><code>def Eulerintegrate(f, y0, tspan):
y = np.zeros([len(tspan),len(y0)])
y[0,:]=y0
for k in range(1, len(tspan)):
y[k,:] = Eulerstep(f, y[k-1], tspan[k]-tspan[k-1])
return y
def RK4integrate(f, y0, tspan):
y = np.zeros([len(tspan),len(y0)])
y[0,:]=y0
for k in range(1, len(tspan)):
y[k,:] = RK4step(f, y[k-1], tspan[k]-tspan[k-1])
return y
</code></pre>
<p>用你给定的问题来调用它们</p>
<pre><code>dt = .1
t = np.arange(0,10,dt)
y0 = np.array([10, 0.0, 10, 10])
sol_euler = Eulerintegrate(orbitsys, y0, t)
x,y,vx,vy = sol_euler.T
plt.plot(x,y)
sol_RK4 = RK4integrate(orbitsys, y0, t)
x,y,vx,vy = sol_RK4.T
plt.plot(x,y)
</code></pre>