Ранг Кутта для системы упругости c маятник

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>

но я не совсем уверен с методом Рутге Кутты. Я сделал симуляцию точки, и это не правильный ответ, что я делаю не так? я думаю, что в оценке 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 не только в качестве постоянной длины покоя пружины, но и в качестве динамической c длины пружины с непредсказуемыми результатами.

В более систематический c подход, можно было бы закодировать систему как систему и использовать векторную версию RK4

def federpendel(u,m,g,l0,k):
    th, r, Vth, Vr = u
    Ath = (-g*np.sin(th) - 2*Vth*Vr) / r
    Ar  = r*Vth**2 - k/m*(r-l0) + g*np.cos(th)
    return np.array([ Vth, Vr, Ath, Ar])

l0 = 10             #spring at rest
g  = 9.81          #gravity
m  = 1             #mass of particle
k  = 40            #spring constant
th0 = 3*np.pi/4 #initial theta
Dth0    = 0         #initial theta velocity
Dr0    = 0         #initial  l velocity
u = np.array([ th0, l0, Dth0, Dr0])
dt = 0.1           #upgrade 
tmax= 20
t = np.arange(0, tmax+0.5*dt, dt)
U = [u];
for n in range(len(t)-1):
    k1 = federpendel(u,m,g,l0,k)*dt
    k2 = federpendel(u+0.5*k1,m,g,l0,k)*dt
    k3 = federpendel(u+0.5*k2,m,g,l0,k)*dt
    k4 = federpendel(u+k3,m,g,l0,k)*dt
    u = u + (k1+2*k2+2*k3+k4)/6
    U.append(u)

th, r, Dth, Dr = np.asarray(U).T
plt.subplot(2,1,1);plt.plot(t,r);
plt.subplot(2,1,2);plt.plot(t,th);
plt.show()

, дающую график