弹性摆系统的runge kutta

Question

我正在尝试使用python解决系统：

<math xmlns="http://www.w3.org/1998/Math/MathML" display="block">
  <mtable columnalign="right left right left right left right left right left right left" rowspacing="3pt" columnspacing="0em 2em 0em 2em 0em 2em 0em 2em 0em 2em 0em" displaystyle="true">
    <mtr>
      <mtd>
        <msub>
          <mrow class="MJX-TeXAtom-ORD">
            <mover>
              <mi>z</mi>
              <mo>&#x02D9;<!-- ˙ --></mo>
            </mover>
          </mrow>
          <mn>1</mn>
        </msub>
      </mtd>
      <mtd>
        <mi></mi>
        <mo>=</mo>
        <mo>&#x2212;<!-- − --></mo>
        <mfrac>
          <mn>1</mn>
          <mi>l</mi>
        </mfrac>
        <mo stretchy="false">[</mo>
        <mi>g</mi>
        <mi>sin</mi>
        <mo>&#x2061;<!-- ⁡ --></mo>
        <mi>&#x03B8;<!-- θ --></mi>
        <mo>+</mo>
        <mn>2</mn>
        <msub>
          <mi>z</mi>
          <mn>1</mn>
        </msub>
        <msub>
          <mi>z</mi>
          <mn>2</mn>
        </msub>
        <mo stretchy="false">]</mo>
        <mo>,</mo>
      </mtd>
    </mtr>
    <mtr>
      <mtd>
        <msub>
          <mrow class="MJX-TeXAtom-ORD">
            <mover>
              <mi>z</mi>
              <mo>&#x02D9;<!-- ˙ --></mo>
            </mover>
          </mrow>
          <mn>2</mn>
        </msub>
      </mtd>
      <mtd>
        <mi></mi>
        <mo>=</mo>
        <mfrac>
          <mn>1</mn>
          <mi>m</mi>
        </mfrac>
        <mo stretchy="false">[</mo>
        <mi>m</mi>
        <mi>l</mi>
        <msubsup>
          <mi>z</mi>
          <mn>1</mn>
          <mn>2</mn>
        </msubsup>
        <mo>&#x2212;<!-- − --></mo>
        <mi>k</mi>
        <mo stretchy="false">(</mo>
        <mi>l</mi>
        <mo>&#x2212;<!-- − --></mo>
        <msub>
          <mi>l</mi>
          <mn>0</mn>
        </msub>
        <mo stretchy="false">)</mo>
        <mo>+</mo>
        <mi>m</mi>
        <mi>g</mi>
        <mi>cos</mi>
        <mo>&#x2061;<!-- ⁡ --></mo>
        <mi>&#x03B8;<!-- θ --></mi>
        <mo stretchy="false">]</mo>
        <mo>.</mo>
      </mtd>
    </mtr>
  </mtable>
</math>

但是我不太确定runge kutta方法。我对此进行了模拟，这不是正确的答案，我在做什么错？我认为ki和mi的评估中存在一些错误，但我读了一百遍，我找不到错误。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle


l0 = 10             #spring at rest
g  = 9.81          #gravity
m  = 1             #mass of particle
k  = 40            #spring constant
dt = 0.1           #upgrade 
Theta0 = 3*np.pi/4 #initial theta
z10    = 0         #initial theta velocity
z20    = 0         #initial  l velocity
tmax, dt = 20, 0.01
t = np.arange(0, tmax+dt, dt)

def f_theta(z1, z2, theta, g, L):
    return (-g*np.sin(theta) - 2*z1*z2) / L

def f_L(z1,theta, g, L, l0, m, k):
    return (m*L*z1**2 - k*(L-l0) + m*g*np.cos(theta)) / m

Thetapoints = []
z1 = []
Lpoints = []
z2 = []

for x in t:
    Thetapoints.append(Theta0)
    z1.append(z10)
    Lpoints.append(l0)
    z2.append(z20)

    m1 = dt*z10
    M1 = dt*f_theta(z10,z20,Theta0,g,l0)

    k1 = dt*z20
    K1 = dt*f_L(z10,Theta0,g,l0,l0,m,k)

    m2 = dt*(z10+0.5*M1)
    M2 = dt*(f_theta(z10+0.5*M1,z20+0.5*K1,Theta0+0.5*m1,g,l0+0.5*k1))

    k2 = dt*(z20+0.5*K1)
    K2 = dt*(f_L(z10+0.5*M2,Theta0+0.5*m2,g,l0+0.5*k2,l0,m,k))

    m3 = dt*(z10+0.5*M2)
    M3 = dt*f_theta(z10+0.5*M2,z20+0.5*K2,Theta0+0.5*m2,g,l0+0.5*k2)

    k3 = dt*(z20+0.5*K2)
    K3 = dt*(f_L(z10+0.5*M2,Theta0+0.5*m2,g,l0+0.5*k2,l0,m,k))

    m4 = dt*(z10+M3)
    M4 = dt*f_theta(z10+M3,z20+K3,Theta0+m3,g,l0+k3)

    k4 = dt*(z20+K3)
    K4 = dt*(f_L(z10+M3,Theta0+m3,g,l0+k3,l0,m,k))

    Theta0 += (m1 + 2*m2 + 2*m3 + m4)/6
    l0     += (k1 + 2*k2 + 2*k3 + k4)/6
    z10    += (M1 + 2*M2 + 2*M3 + M4)/6
    z20    += (K1 + 2*K2 + 2*K3 + K4)/6

x =  np.array(Lpoints)* np.sin(np.array(Thetapoints))
y = -np.array(Lpoints)* np.cos(np.array(Thetapoints))

plt.plot(t,Lpoints)
plt.show()

Answer 1

最明显的错误是，您不仅将l0用作弹簧的恒定静止长度，而且还将其用作弹簧的动态长度，其结果无法预测。