runge kutta for system of elastic pendulum

Question

I'm trying to solve the system with 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>

but i'm not quite sure with the runge kutta method. I made a simulation of the point and this is not the right answer, what i am i doing wrong? i think there is some mistake in the evaluation of ki and mi but i read it hundred times and i can't find the mistake.

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()

You should add that z1 is the velocity of theta and z2 the velocity of the length. Your naming scheme and parameter list is a little confusing. — Lutz Lehmann, Jan 16 '20 at 08:40
Yes, with some additional context, the equations are correct. What is wrong is the use of `l0`. The rest length should be a constant, not at the same time the length variable. If you decouple this use (and the step size is small enough), you should get the correct results. — Lutz Lehmann, Jan 16 '20 at 08:44

score 2 · Answer 1 · answered Jan 16 '20 at 09:49

The most obvious error is that you use l0 not only as the constant rest length of the spring, but also as the dynamic length of the spring, with unpredictable results.

In a more systematic approach, one would encode the system as system and use the vector version of 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()

giving the plot

Very nice approach, Surely more compact and immediate. I will implement it but for the first time i wanted to have a little more "control" over the parameter. But surely i will follow the advice! — Salvatore Lampitelli, Jan 16 '20 at 18:35

runge kutta for system of elastic pendulum

1 Answers1