Pendulums¶

CalcpadCE handles pendulum dynamics end-to-end on a single worksheet: nonlinear ODEs are integrated numerically, while inline SVG turns each time step into a frame of a physics animation.

All five sheets share the same recipe — derive the equations of motion, integrate with Runge-Kutta 4, then draw the bob, rod and trajectory inside the worksheet. The simple undamped and simple damped cases stay in one degree of freedom, so the integration result can be compared against analytical and small-angle approximations on the same plot — see the side-by-side comparison sheet for that. The elastic-damped variant adds a radial spring degree of freedom, and the double pendulum couples two angles into the chaotic regime.

Lengths, masses, damping factors and initial conditions are exposed at the top of each sheet, so every animation doubles as a parametric playground for stability and energy-loss studies.

Simple Undamped Pendulum 🎬¶

Single-degree-of-freedom undamped pendulum integrated with classical RK4 from a near-inverted initial angle. Each time step is rendered as an SVG frame so the swing builds up live, alongside an angle-versus-time plot.

Code:

#rad
'<table><tr><td>
'<h4>Input parameters</h4>
'Gravitational acceleration (m/s²) - 'g = 9.80665m/s^2
'Pendulum length -'l = 0.5m
'Pendulum mass -'m = 1kg
'Initial angle -'θ_0 = 179°|rad
'Maximum simulation time -'t_max = 11s
'<h4>Theoretical solution for large rotations</h4>
'</td><td>
#hide
Precision = 10^-5
PlotWidth = 300','PlotHeight = 150','PlotSVG = 1
w = 1.2*l', 'h_ = 2/3*w
W = 150', 'H = 130
k = W/w
L = l*k
R = L/10
x = -L*sin(θ_0)
y = L*cos(θ_0)
δ = 16
δ_x = 0.25*δ
δ_y = -2*δ
y_g = y - δ_y + R
#def svg$ = '<svg viewbox="'-60 - W/2' '-H-20' 'W + 120' '2*H + 40'" xmlns="http://www.w3.org/2000/svg" version="1.1" style="font-family: Georgia Pro; font-size:14px; width:'W'pt; height:'H + 50'pt; margin-left:20pt;">
#def thin_style$ = style="stroke:goldenrod; stroke-width:1; fill:none"
#def thick_style$ = style="stroke:steelblue; stroke-width:2; fill:none"
#def ball_style$ = style="stroke:steelblue; stroke-width:0.5;  fill:url(#ball)"
#def solid_style$ = style="stroke:black; stroke-width:1; fill:#eee"
#def load_style$ = style="stroke:steelblue; stroke-width:1; fill:steelblue;"
#include svg_drawing.cpd
#show
#val
svg$
'<defs><radialGradient id="ball" cx="35%" cy="35%"><stop offset="0%" stop-color="lightcyan"/><stop offset="100%" stop-color="lightblue"/></radialGradient><pattern id="concrete" x="4" y="4" width="8" height="8" patternUnits="userSpaceOnUse"><rect x="0" y="0" width="8" height="8" fill="#eae9e8" /><circle cx="1" cy="1" r="1.2" fill="#ccb" /><circle cx="5" cy="2" r="1.6" fill="#dadad0" /><circle cx="4" cy="6" r="0.8" fill="#aa9" /><circle cx="3" cy="4" r="0.4" fill="#884" /><circle cx="7" cy="5" r="1.2" fill="#cacaba" /><circle cx="5" cy="3" r="0.9" fill="#fffded" /></pattern></defs>
'<polygon points="0,0 '-δ','-δ' 'δ','-δ'" solid_style$/>
line$(0; 0; 0; l*k; thin_style$)
'<path d="M 0 '0.2*L' A '0.2*L' '0.2*L' 0 0 1 '0.2*x' '0.2*y'" thin_style$/>
texth$(δ*sign(θ_0);0.35*L;'θ_0*180/π'°)
line$(0; 0; x; y; thick_style$)
text$(x/2-δ_x*y/L; y/2+δ_x*x/L;450°-θ_0;'l'm)
circle$(0; 0; R/2; solid_style$)
circle$(0; 0; R/8; solid_style$)
circle$(x; y; R; ball_style$)
'<polygon points="'x','y_g + δ_y' 'x','y_g' 'x - δ_x','y_g + δ_y/3' 'x + δ_x','y_g + δ_y/3' 'x','y_g'" load_style$/>
texth$(x+3*δ_x; y_g-δ;g)
'</svg>
#equ
'</td></tr></table>
#noc
'Differential equation -'θ″ + g/l*sin(θ) = 0
#equ
'Incomplete elliptic integral of the first kind
F(φ; k) = $Integral{1/sqrt(1 - k^2*sin(θ)^2) @ θ = 0 : φ}
'Jacobi elliptic functions
'Modulus -'k = sin(θ_0/2)
am(u; k) = $Root{F(φ; k) = u @ φ = 0 : 3*π}
sn(u; k) = sin(am(u; k))','cn(u; k) = cos(am(u; k))
dn(u; k) = sqrt(1 - k*sn(u; k)^2)','cd(u; k) = cn(u; k)/dn(u; k)
'Period -'T_e = 4*sqrt(l/g)*F(π/2; k)
'Cyclic frequency -'f_e = 1/T_e|Hz
'Angular frequency -'ω_e = 2*π*rad*f_e|rad/s
'Equation of motion -'θ_e(t) = 2*asin(k*cd(sqrt(g/l)*t; k))
'Energy -'E_0 = m*l*g*(1 - cos(θ_0))|J
'<h4>Solution by Runge-Kutta RK4 method (explicit)</h4>
'For that purpose, the II order equation is reduced to the following system of I order equations
#noc
θ′ = ω'and'ω′ = (g/l)*sin(-θ)
#equ
'The solution will be performed iteratively
'Step size -'h = 0.02s
'Number of steps -'n = t_max/h
#hide
N = min(n; 80)
k_N = floor(n/N)
N = n\k_N
#show
'For each time step'n'the values for the next step will be obtained by using the following equstions
#noc
'<table class="eq" style="margin-left:24pt;">
'<tr><td>First step (k₁) -</td><td>'k_1,θ = ω_i', </td><td>' _
k_1,ω = (g/l)*sin(-θ_i)'</td></tr>
'<tr><td>Second step (k₂) -</td><td>'k_2,θ = ω_i + 0.5*h*k_1,ω', </td><td>'k_2,ω = (g/l)*sin(-(θ_i + 0.5*h*k_1,θ))'</td></tr>
'<tr><td>Third step (k₃) -</td><td>'k_3,θ = ω_i + 0.5*h*k_2,ω', </td><td>'k_3,ω = (g/l)*sin(-(θ_i + 0.5*h*k_2,θ))'</td></tr>
'<tr><td>Fourth step (k₄) -</td><td>'k_4,θ = ω_i + h*k_3,ω', </td><td>'k_4,ω = (g/l)*sin(-(θ_i + h*k_3,θ))'</td></tr>
'</table>
'Update values using weighted averages
'<p class="eq" style="border-left:solid 1pt black; margin-left:24pt; padding-left:3pt;"><var>θ</var><sub>n+1</sub> ='θ_i + (h/6)*(k_1,θ + 2*k_2,θ + 2*k_3,θ + k_4,θ)'<br/> <var>ω</var><sub>n+1</sub> = 'ω_i + (h/6)*(k_1,ω + 2*k_2,ω + 2*k_3,ω + k_4,ω)'</p>
#equ
'Allocate vectors
θ_RK4 = vector(n)'<br/>'ω_RK4 = vector(n)'<br/>'E_RK4 = vector(n)
'Set initial conditions
θ_RK4.1 = θ_0/1rad','ω_RK4.1 = 0/s
'Perform Runge-Kutta 4 steps
#for i = 1 : n - 1
    'RK4 factors
    k_1,θ = ω_RK4.i', ' _
    k_1,ω = (g/l)*sin(-θ_RK4.i)
    k_2,θ = ω_RK4.i + 0.5*h*k_1,ω
    k_2,ω = (g/l)*sin(-(θ_RK4.i + 0.5*h*k_1,θ))
    k_3,θ = ω_RK4.i + 0.5*h*k_2,ω
    k_3,ω = (g/l)*sin(-(θ_RK4.i + 0.5*h*k_2,θ))
    k_4,θ = ω_RK4.i + h*k_3,ω
    k_4,ω = (g/l)*sin(-(θ_RK4.i + h*k_3,θ))
    'Update values using weighted averages
    'Rotation -'θ_RK4.(i + 1) = θ_RK4.i + (h/6)*(k_1,θ + 2*k_2,θ + 2*k_3,θ + k_4,θ)
    'Angular velocity -'ω_RK4.(i + 1) = ω_RK4.i + (h/6)*(k_1,ω + 2*k_2,ω + 2*k_3,ω + k_4,ω)
    'Energy -'E_RK4.i = m*l^2*(1/2*ω_RK4.i^2 + (g/l)*(1 - cos(θ_RK4.i)))|J
    #hide
#loop
E_RK4.n = m*l^2*(1/2*ω_RK4.n^2 + (g/l)*(1 - cos(θ_RK4.n)))|J
t(i) = i*h - h
r2d = 180/π
#show
'. . .
'<style>.Series2{stroke:Lime!Important; stroke-dasharray: 30, 10;}</style>
'<h4>Plot results</h4>
'<table><tr><td>
'Rotation - <var>θ</var>(<var>t</var>), [°]<br/>
$Plot{t(i)|θ_e(t(i))*r2d & t(i)|θ_RK4.i*r2d @ i = 1 : min(1.5*N; n)}
'</td><td><br/><br/><br/>
'<b style="color:Red">━━━</b> Theoretical <br/>
'<b style="color:Lime">╸╸╸</b> Runge-Kutta 4
'</td></tr><tr><td>
'Absolute error Δθ(<var>t</var>), [°]<br/>
$Plot{t(i)|(θ_RK4.i - θ_e(t(i)))*r2d @ i = 1 : min(1.5*N; n)}
'</td><td><br/><br/><br/>
'<p><b style="color:Red">━━━</b> Runge-Kutta 4
'</td></tr><tr><td>
'Energy <var>E</var>(<var>t</var>), [J]<br/>
$Plot{t(i)|E_0 & t(i)|E_RK4.i & 0s|0.99999*E_0 @ i = 1 : min(1.5*N; n)}
'</td><td><br/><br/><br/>
'<b style="color:Red">━━━</b> Theoretical <br/>
'<b style="color:Lime">╸╸╸</b> Runge-Kutta 4
'</td></tr></table>
'<style>.fr{display:none;} .fr svg{margin-right:64pt; align:center;}</style>
'<h4>Animate results</h4>
#format 0.00
#hide
k = W/w
#val
#for i_N = 1 : N
    #hide
    i = i_N*k_N
    θ = θ_RK4.i
    x = -L*sin(θ)
    y = L*cos(θ)
    t = t(i)
    #show
    #val
    '<div class="fr" id="p1-fr-'i_N'">
    #equ
    #nosub
    'Simulation time -'t'<br/>
    $Plot{t(j)|θ_RK4.j*r2d & t(i)|θ_RK4.i*r2d & t_max|θ_0/1rad*r2d & t_max|-θ_0/1rad*r2d @ j = 1 : i}
    #val
    svg$
    '<polygon points="0,0 '-δ','-δ' 'δ','-δ'" solid_style$/>
    line$(0; 0; x; y; thick_style$)
    circle$(0; 0; R/2; solid_style$)
    circle$(0; 0; R/8; solid_style$)
    circle$(x; y; R; ball_style$)
    '</svg>
    '</div>
#loop
'<script>(function(){$("#p1-fr-1").show();var i=1;var fr=$("#p1-fr-1");setInterval(function(){fr.hide();if(++i>'N')i=1;fr=$("#p1-fr-"+i);fr.show();},'1000*h*k_N');})();</script>
#equ

Rendered Output:

Input parameters

Gravitational acceleration (m/s²) - g = 9.81 ms²

Pendulum length - l = 0.5 m

Pendulum mass - m = 1 kg

Initial angle - θ₀ = 179° = 3.12 rad

Maximum simulation time - t_max = 11 s

Theoretical solution for large rotations

Differential equation - θ″ + gl · sin ( θ ) = 0

Incomplete elliptic integral of the first kind

F ( φ; k ) = φ∫0 1 √ 1 − k² · sin ( θ ) ² dθ

Jacobi elliptic functions

Modulus - k = sin(θ₀2) = sin(3.12 rad2) = 1

am ( u; k ) = $Root{F ( φ; k ) = u; φ ∈ [0; 3 · π]}

sn ( u; k ) = sin ( am ( u; k ) ) , cn ( u; k ) = cos ( am ( u; k ) )

dn ( u; k ) = √ 1 − k · sn ( u; k ) ² , cd ( u; k ) = cn ( u; k ) dn ( u; k )

Period - T_e = 4 · lg · F (π2; k) = 4 · 0.5 m9.81 m ∕ s² · F (3.142; 1) = 5.53 s

Cyclic frequency - f_e = 1T_e = 15.53 s = 0.181 Hz

Angular frequency - ω_e = 2 · π · rad · f_e = 2 · 3.14 · rad · 0.181 Hz = 1.14 rad ∕ s

Equation of motion - θ_e ( t ) = 2 · asin(k · cd ( gl · t; k))

Energy - E₀ = m · l · g · ( 1 − cos ( θ₀ ) ) = 1 kg · 0.5 m · 9.81 m ∕ s² · ( 1 − cos ( 3.12 rad ) ) = 9.81 J

Solution by Runge-Kutta RK4 method (explicit)

For that purpose, the II order equation is reduced to the following system of I order equations

θ′ = ω and ω′ = gl · sin ( -θ )

The solution will be performed iteratively

Step size - h = 0.02 s

Number of steps - n = t_maxh = 11 s0.02 s = 550

For each time step n = 550 the values for the next step will be obtained by using the following equstions

First step (k₁) -	`k`_1,θ = `ω`_i ,	`k`_1,ω = `gl` · sin ( -`θ`_i )
Second step (k₂) -	`k`_2,θ = `ω`_i + 0.5 · `h` · `k`_1,ω ,	`k`_2,ω = `gl` · sin ( - ( `θ`_i + 0.5 · `h` · `k`_1,θ ) )
Third step (k₃) -	`k`_3,θ = `ω`_i + 0.5 · `h` · `k`_2,ω ,	`k`_3,ω = `gl` · sin ( - ( `θ`_i + 0.5 · `h` · `k`_2,θ ) )
Fourth step (k₄) -	`k`_4,θ = `ω`_i + `h` · `k`_3,ω ,	`k`_4,ω = `gl` · sin ( - ( `θ`_i + `h` · `k`_3,θ ) )

Update values using weighted averages

θ_n+1 = θ_i + h6 · ( k_1,θ + 2 · k_2,θ + 2 · k_3,θ + k_4,θ )
ω_n+1 = ω_i + h6 · ( k_1,ω + 2 · k_2,ω + 2 · k_3,ω + k_4,ω )

Allocate vectors

⃗θ_RK4 = vector ( n ) = vector ( 550 ) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... 0]
⃗ω_RK4 = vector ( n ) = vector ( 550 ) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... 0]
⃗E_RK4 = vector ( n ) = vector ( 550 ) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... 0]

Set initial conditions

⃗θ_RK4.1 = θ₀1 rad = 3.12 rad1 rad = 3.12 , ⃗ω_RK4.1 = 0s = 0 s^-1

Perform Runge-Kutta 4 steps

RK4 factors

k_1,θ = ⃗ω_RK4.1 = 0 s^-1 , k_1,ω = gl · sin ( -⃗θ_RK4.1 ) = 9.81 m ∕ s²0.5 m · sin ( -3.12 ) = -0.342 s^-2

k_2,θ = ⃗ω_RK4.1 + 0.5 · h · k_1,ω = 0 s^-1 + 0.5 · 0.02 s · ( -0.342 s^-2 ) = -0.00342 s^-1

k_2,ω = gl · sin ( - ( ⃗θ_RK4.1 + 0.5 · h · k_1,θ ) ) = 9.81 m ∕ s²0.5 m · sin ( - ( 3.12 + 0.5 · 0.02 s · 0 s^-1 ) ) = -0.342 s^-2

k_3,θ = ⃗ω_RK4.1 + 0.5 · h · k_2,ω = 0 s^-1 + 0.5 · 0.02 s · ( -0.342 s^-2 ) = -0.00342 s^-1

k_3,ω = gl · sin ( - ( ⃗θ_RK4.1 + 0.5 · h · k_2,θ ) ) = 9.81 m ∕ s²0.5 m · sin ( - ( 3.12 + 0.5 · 0.02 s · ( -0.00342 s^-1 ) ) ) = -0.343 s^-2

k_4,θ = ⃗ω_RK4.1 + h · k_3,ω = 0 s^-1 + 0.02 s · ( -0.343 s^-2 ) = -0.00686 s^-1

k_4,ω = gl · sin ( - ( ⃗θ_RK4.1 + h · k_3,θ ) ) = 9.81 m ∕ s²0.5 m · sin ( - ( 3.12 + 0.02 s · ( -0.00342 s^-1 ) ) ) = -0.344 s^-2

Update values using weighted averages

Rotation - ⃗θ_RK4.2 = ⃗θ_RK4.1 + h6 · ( k_1,θ + 2 · k_2,θ + 2 · k_3,θ + k_4,θ ) = 3.12 + 0.02 s6 · ( 0 s^-1 + 2 · ( -0.00342 s^-1 ) + 2 · ( -0.00342 s^-1 ) + ( -0.00686 s^-1 ) ) = 3.12

Angular velocity - ⃗ω_RK4.2 = ⃗ω_RK4.1 + h6 · ( k_1,ω + 2 · k_2,ω + 2 · k_3,ω + k_4,ω ) = 0 s^-1 + 0.02 s6 · ( -0.342 s^-2 + 2 · ( -0.342 s^-2 ) + 2 · ( -0.343 s^-2 ) + ( -0.344 s^-2 ) ) = -0.00685 s^-1

Energy - ⃗E_RK4.1 = m · l² · (12 · ⃗ω_RK4.1² + gl · ( 1 − cos ( ⃗θ_RK4.1 ) ) ) = 1 kg · ( 0.5 m ) ² · (12 · 0 s^-1² + 9.81 m ∕ s²0.5 m · ( 1 − cos ( 3.12 ) ) ) = 9.81 J

. . .

Plot results

Rotation - `θ`(`t`), [°]	━━━ Theoretical ╸╸╸ Runge-Kutta 4
Absolute error Δθ(`t`), [°]	━━━ Runge-Kutta 4
Energy `E`(`t`), [J]	━━━ Theoretical ╸╸╸ Runge-Kutta 4

Animate results

Simulation time - t = 0.10 s

Simulation time - t = 0.22 s

Simulation time - t = 0.34 s

Simulation time - t = 0.46 s

Simulation time - t = 0.58 s

Simulation time - t = 0.70 s

Simulation time - t = 0.82 s

Simulation time - t = 0.94 s

Simulation time - t = 1.06 s

Simulation time - t = 1.18 s

Simulation time - t = 1.30 s

Simulation time - t = 1.42 s

Simulation time - t = 1.54 s

Simulation time - t = 1.66 s

Simulation time - t = 1.78 s

Simulation time - t = 1.90 s

Simulation time - t = 2.02 s

Simulation time - t = 2.14 s

Simulation time - t = 2.26 s

Simulation time - t = 2.38 s

Simulation time - t = 2.50 s

Simulation time - t = 2.62 s

Simulation time - t = 2.74 s

Simulation time - t = 2.86 s

Simulation time - t = 2.98 s

Simulation time - t = 3.10 s

Simulation time - t = 3.22 s

Simulation time - t = 3.34 s

Simulation time - t = 3.46 s

Simulation time - t = 3.58 s

Simulation time - t = 3.70 s

Simulation time - t = 3.82 s

Simulation time - t = 3.94 s

Simulation time - t = 4.06 s

Simulation time - t = 4.18 s

Simulation time - t = 4.30 s

Simulation time - t = 4.42 s

Simulation time - t = 4.54 s

Simulation time - t = 4.66 s

Simulation time - t = 4.78 s

Simulation time - t = 4.90 s

Simulation time - t = 5.02 s

Simulation time - t = 5.14 s

Simulation time - t = 5.26 s

Simulation time - t = 5.38 s

Simulation time - t = 5.50 s

Simulation time - t = 5.62 s

Simulation time - t = 5.74 s

Simulation time - t = 5.86 s

Simulation time - t = 5.98 s

Simulation time - t = 6.10 s

Simulation time - t = 6.22 s

Simulation time - t = 6.34 s

Simulation time - t = 6.46 s

Simulation time - t = 6.58 s

Simulation time - t = 6.70 s

Simulation time - t = 6.82 s

Simulation time - t = 6.94 s

Simulation time - t = 7.06 s

Simulation time - t = 7.18 s

Simulation time - t = 7.30 s

Simulation time - t = 7.42 s

Simulation time - t = 7.54 s

Simulation time - t = 7.66 s

Simulation time - t = 7.78 s

Simulation time - t = 7.90 s

Simulation time - t = 8.02 s

Simulation time - t = 8.14 s

Simulation time - t = 8.26 s

Simulation time - t = 8.38 s

Simulation time - t = 8.50 s

Simulation time - t = 8.62 s

Simulation time - t = 8.74 s

Simulation time - t = 8.86 s

Simulation time - t = 8.98 s

Simulation time - t = 9.10 s

Simulation time - t = 9.22 s

Simulation time - t = 9.34 s

Simulation time - t = 9.46 s

Simulation time - t = 9.58 s

Simulation time - t = 9.70 s

Simulation time - t = 9.82 s

Simulation time - t = 9.94 s

Simulation time - t = 10.06 s

Simulation time - t = 10.18 s

Simulation time - t = 10.30 s

Simulation time - t = 10.42 s

Simulation time - t = 10.54 s

Simulation time - t = 10.66 s

Simulation time - t = 10.78 s

Simulation time - t = 10.90 s

Simple Damped Pendulum 🎬¶

Adds linear viscous damping $b = β\, b_{cr}$ to the simple pendulum and integrates the resulting nonlinear ODE with RK4. The animation shows amplitude decay toward equilibrium for an arbitrary damping ratio.

Code:

#rad
'<table><tr><td>
'<h4>Input parameters</h4>
'Gravitational acceleration -'g = 9.80665m/s^2
'Pendulum length -'l = 1m
'Pendulum mass -'m = 1kg
'Initial angle -'θ_0 = -179°|rad
'Initial angular velocity -'ω_0 = 1.834rad/s
'Maximum simulation time -'t_max = 20s
'Damping factor (normaized) -'β = 0.005
'Critical damping -'b_cr = 2*m*sqrt(g/l)
'Damping -'b = β*b_cr
'</td><td>
#hide
Precision = 10^-5
PlotWidth = 300','PlotHeight = 150','PlotSVG = 1
w = 1.2*l', 'h_ = 2/3*w
W = 180', 'H = 150
k = W/w
L = l*k
R = L/10
x = L*sin(θ_0)
y = L*cos(θ_0)
δ = 16
δ_x = 0.25*δ
δ_y = -2*δ
y_g = y - δ_y + R
#def svg$ = '<svg viewbox="'-20-W/2' '-H-20' 'W+40' '2*H+40'" xmlns="http://www.w3.org/2000/svg" version="1.1" style="font-family: Georgia Pro; font-size:14px; width:'W+50'pt; height:'H+50'pt; margin-left:20pt;">
#def thin_style$ = style="stroke:goldenrod; stroke-width:1; fill:none"
#def thick_style$ = style="stroke:steelblue; stroke-width:2; fill:none"
#def ball_style$ = style="stroke:steelblue; stroke-width:0.5;  fill:url(#ball)"
#def solid_style$ = style="stroke:black; stroke-width:1; fill:#eee"
#def load_style$ = style="stroke:steelblue; stroke-width:1; fill:steelblue;"
#include svg_drawing.cpd
#show
#val
svg$
'<defs><radialGradient id="ball" cx="35%" cy="35%"><stop offset="0%" stop-color="lightcyan"/><stop offset="100%" stop-color="lightblue"/></radialGradient><pattern id="concrete" x="4" y="4" width="8" height="8" patternUnits="userSpaceOnUse"><rect x="0" y="0" width="8" height="8" fill="#eae9e8" /><circle cx="1" cy="1" r="1.2" fill="#ccb" /><circle cx="5" cy="2" r="1.6" fill="#dadad0" /><circle cx="4" cy="6" r="0.8" fill="#aa9" /><circle cx="3" cy="4" r="0.4" fill="#884" /><circle cx="7" cy="5" r="1.2" fill="#cacaba" /><circle cx="5" cy="3" r="0.9" fill="#fffded" /></pattern></defs>
'<polygon points="0,0 '-δ','-δ' 'δ','-δ'" solid_style$/>
line$(0; 0; 0; l*k; thin_style$)
'<path d="M 0 '0.2*L' A '0.2*L' '0.2*L' 0 0 'θ_0 < 0°' '0.2*x' '0.2*y'" thin_style$/>
texth$(-1.2*δ*sign(θ_0);0.25*L;'θ_0*180/π'°)
line$(0; 0; x; y; thick_style$)
text$(x/2-δ_x*y/L; y/2+δ_x*x/L;if(θ_0<0°;450°+θ_0;450°-θ_0);'l'm)
circle$(0; 0; R/2; solid_style$)
circle$(0; 0; R/8; solid_style$)
circle$(x; y; R; ball_style$)
'<polygon points="'x','y_g + δ_y' 'x','y_g' 'x - δ_x','y_g + δ_y/3' 'x + δ_x','y_g + δ_y/3' 'x','y_g'" load_style$/>
texth$(x+3*δ_x; y_g-δ;g)
'</svg>
#equ
'</td></tr></table>
'<h4>Solution by Runge-Kutta RK4 method (explicit)</h4>
#noc
'Differential equation -'θ″ + b/m*θ′ + g/l*sin(θ) = 0
#equ
'The II order equation is first reduced to the following system of I order equations
#noc
θ′ = ω'and'ω′ = (g/l)*sin(-θ) - b/m*ω
#equ
'The solution is performed iteratively
'Step size -'h = 0.02s
'Number of steps -'n = t_max/h
#hide
N = min(n; 120)
k_N = floor(n/N)
N = n\k_N
#show
'For each time step'n'the values for the next step will be obtained by using the following equstions
#noc
'<table class="eq" style="margin-left:24pt;">
'<tr><td>First step (k₁) -</td><td>'k_1,θ = ω_i', </td><td>' _
k_1,ω = (g/l)*sin(-θ_i) - b/m*k_1,θ'</td></tr>
'<tr><td>Second step (k₂) -</td><td>'k_2,θ = ω_i + 0.5*h*k_1,ω', </td><td>'k_2,ω = (g/l)*sin(-(θ_i + 0.5*h*k_1,θ)) - b/m*k_2,θ'</td></tr>
'<tr><td>Third step (k₃) -</td><td>'k_3,θ = ω_i + 0.5*h*k_2,ω', </td><td>'k_3,ω = (g/l)*sin(-(θ_i + 0.5*h*k_2,θ)) - b/m*k_3,θ'</td></tr>
'<tr><td>Fourth step (k₄) -</td><td>'k_4,θ = ω_i + h*k_3,ω', </td><td>'k_4,ω = (g/l)*sin(-(θ_i + h*k_3,θ)) - b/m*k_4,θ'</td></tr>
'</table>
'Update values using weighted averages
'<p class="eq" style="border-left:solid 1pt black; margin-left:24pt; padding-left:3pt;"><var>θ</var><sub>n+1</sub> ='θ_i + (h/6)*(k_1,θ + 2*k_2,θ + 2*k_3,θ + k_4,θ)'<br/> <var>ω</var><sub>n+1</sub> = 'ω_i + (h/6)*(k_1,ω + 2*k_2,ω + 2*k_3,ω + k_4,ω)'</p>
#equ
'Allocate vectors
θ = vector(n)'<br/>'ω = vector(n)
'Set initial conditions
θ.1 = θ_0/1rad','ω.1 = ω_0/1rad
'Perform Runge-Kutta 4 steps
#for i = 1 : n - 1
    'RK4 factors
    k_1,θ = ω.i', ' _
    k_1,ω = (g/l)*sin(-θ.i) - b/m*k_1,θ
    k_2,θ = ω.i + 0.5*h*k_1,ω
    k_2,ω = (g/l)*sin(-(θ.i + 0.5*h*k_1,θ)) - b/m*k_2,θ
    k_3,θ = ω.i + 0.5*h*k_2,ω
    k_3,ω = (g/l)*sin(-(θ.i + 0.5*h*k_2,θ)) - b/m*k_3,θ
    k_4,θ = ω.i + h*k_3,ω
    k_4,ω = (g/l)*sin(-(θ.i + h*k_3,θ)) - b/m*k_4,θ
    'Update values using weighted averages
    'Rotation -'θ.(i + 1) = θ.i + (h/6)*(k_1,θ + 2*k_2,θ + 2*k_3,θ + k_4,θ)
    'Angular velocity -'ω.(i + 1) = ω.i + (h/6)*(k_1,ω + 2*k_2,ω + 2*k_3,ω + k_4,ω)' < /p >
    #hide
#loop
t(i) = i*h - h
r2d = 180/π
#show
'Kinetic energy -'E_k(i) = m/2*l^2*ω.i^2|J
'Potential energy-'E_p(i) = m*l*g*(1 - cos(θ.i))|J
'...
'<h4>Plot results</h4>
'Rotation - <var>θ</var>(<var>t</var>), [°]
$Plot{t(i)|θ.i*r2d @ i = 1 : n}
'Energy - <var>E</var>(<var>t</var>), [J]
$Plot{t(i)|E_k(i) + E_p(i) & 0s|0J & 0s|0J & 0s|0J & t(i)|E_p(i) & t(i)|E_k(i) & 0s|0J @ i = 1 : n}
'<br/><b style="color:MediumVioletRed">━━</b> Potential - <var>E</var><sub>p</sub>,
'<b style="color:MediumSpringGreen">━━</b> Kinetic - <var>E</var><sub>k</sub>,
'<b style="color:Tomato">━━</b> Total - <var>E</var><sub>tot</sub>
'<style>.Series3, .Series4{stroke:none!Important; fill:none!Important;} .fr{display:none;} .fr svg{margin-right:64pt; align-center;}</style>
'<h4>Animate results</h4>
#format 0.00
#hide
k = W/w
#val
#for i_N = 1 : N
    #hide
    i = i_N*k_N
    θ_i = θ.i
    x = L*sin(θ_i)
    y = L*cos(θ_i)
    t = t(i)
    #show
    #val
    '<div class="fr" id="p2-fr-'i_N'">
    #equ
    #nosub
    'Simulation time -'t'<br/>
    $Plot{t(i)|θ.i*r2d & t(j)|θ.j*r2d & t_max|θ_0/1rad*r2d & t_max|-θ_0/1rad*r2d @ j = 1 : i}
    #val
    svg$
    '<polygon points="0,0 '-δ','-δ' 'δ','-δ'" solid_style$/>
    line$(0; 0; x; y; thick_style$)
    circle$(0; 0; R/2; solid_style$)
    circle$(0; 0; R/8; solid_style$)
    circle$(x; y; R; ball_style$)
    '</svg>
    '</div>
#loop
'<script>(function(){$("#p2-fr-1").show();var i=1;var fr=$("#p2-fr-1");setInterval(function(){fr.hide();if(++i>'N')i=1;fr=$("#p2-fr-"+i);fr.show();},'1000*h*k_N');})();</script>
#equ

Rendered Output:

Input parameters

Gravitational acceleration - g = 9.81 ms²

Pendulum length - l = 1 m

Pendulum mass - m = 1 kg

Initial angle - θ₀ = -179° = -3.12 rad

Initial angular velocity - ω₀ = 1.83 rads

Maximum simulation time - t_max = 20 s

Damping factor (normaized) - β = 0.005

Critical damping - b_cr = 2 · m · gl = 2 · 1 kg · 9.81 m ∕ s²1 m = 6.26 kg ∕ s

Damping - b = β · b_cr = 0.005 · 6.26 kg ∕ s = 0.0313 kg ∕ s

Solution by Runge-Kutta RK4 method (explicit)

Differential equation - θ″ + bm · θ′ + gl · sin ( θ ) = 0

The II order equation is first reduced to the following system of I order equations

θ′ = ω and ω′ = gl · sin ( -θ ) − bm · ω

The solution is performed iteratively

Step size - h = 0.02 s

Number of steps - n = t_maxh = 20 s0.02 s = 1000

For each time step n = 1000 the values for the next step will be obtained by using the following equstions

First step (k₁) -	`k`_1,θ = `ω`_i ,	`k`_1,ω = `gl` · sin ( -`θ`_i ) − `bm` · `k`_1,θ
Second step (k₂) -	`k`_2,θ = `ω`_i + 0.5 · `h` · `k`_1,ω ,	`k`_2,ω = `gl` · sin ( - ( `θ`_i + 0.5 · `h` · `k`_1,θ ) ) − `bm` · `k`_2,θ
Third step (k₃) -	`k`_3,θ = `ω`_i + 0.5 · `h` · `k`_2,ω ,	`k`_3,ω = `gl` · sin ( - ( `θ`_i + 0.5 · `h` · `k`_2,θ ) ) − `bm` · `k`_3,θ
Fourth step (k₄) -	`k`_4,θ = `ω`_i + `h` · `k`_3,ω ,	`k`_4,ω = `gl` · sin ( - ( `θ`_i + `h` · `k`_3,θ ) ) − `bm` · `k`_4,θ

Update values using weighted averages

θ_n+1 = θ_i + h6 · ( k_1,θ + 2 · k_2,θ + 2 · k_3,θ + k_4,θ )
ω_n+1 = ω_i + h6 · ( k_1,ω + 2 · k_2,ω + 2 · k_3,ω + k_4,ω )

Allocate vectors

⃗θ = vector ( n ) = vector ( 1000 ) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... 0]
⃗ω = vector ( n ) = vector ( 1000 ) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... 0]

Set initial conditions

⃗θ₁ = θ₀1 rad = -3.12 rad1 rad = -3.12 , ⃗ω₁ = ω₀1 rad = 1.83 rad ∕ s1 rad = 1.83 s^-1

Perform Runge-Kutta 4 steps

RK4 factors

k_1,θ = ⃗ω₁ = 1.83 s^-1 , k_1,ω = gl · sin ( -⃗θ₁ ) − bm · k_1,θ = 9.81 m ∕ s²1 m · sin ( - ( -3.12 ) ) − 0.0313 kg ∕ s1 kg · 1.83 s^-1 = 0.114 s^-2

k_2,θ = ⃗ω₁ + 0.5 · h · k_1,ω = 1.83 s^-1 + 0.5 · 0.02 s · 0.114 s^-2 = 1.84 s^-1

k_2,ω = gl · sin ( - ( ⃗θ₁ + 0.5 · h · k_1,θ ) ) − bm · k_2,θ = 9.81 m ∕ s²1 m · sin ( - ( -3.12 + 0.5 · 0.02 s · 1.83 s^-1 ) ) − 0.0313 kg ∕ s1 kg · 1.84 s^-1 = 0.293 s^-2

k_3,θ = ⃗ω₁ + 0.5 · h · k_2,ω = 1.83 s^-1 + 0.5 · 0.02 s · 0.293 s^-2 = 1.84 s^-1

k_3,ω = gl · sin ( - ( ⃗θ₁ + 0.5 · h · k_2,θ ) ) − bm · k_3,θ = 9.81 m ∕ s²1 m · sin ( - ( -3.12 + 0.5 · 0.02 s · 1.84 s^-1 ) ) − 0.0313 kg ∕ s1 kg · 1.84 s^-1 = 0.294 s^-2

k_4,θ = ⃗ω₁ + h · k_3,ω = 1.83 s^-1 + 0.02 s · 0.294 s^-2 = 1.84 s^-1

k_4,ω = gl · sin ( - ( ⃗θ₁ + h · k_3,θ ) ) − bm · k_4,θ = 9.81 m ∕ s²1 m · sin ( - ( -3.12 + 0.02 s · 1.84 s^-1 ) ) − 0.0313 kg ∕ s1 kg · 1.84 s^-1 = 0.474 s^-2

Update values using weighted averages

Rotation - ⃗θ₂ = ⃗θ₁ + h6 · ( k_1,θ + 2 · k_2,θ + 2 · k_3,θ + k_4,θ ) = -3.12 + 0.02 s6 · ( 1.83 s^-1 + 2 · 1.84 s^-1 + 2 · 1.84 s^-1 + 1.84 s^-1 ) = -3.09

Angular velocity - ⃗ω₂ = ⃗ω₁ + h6 · ( k_1,ω + 2 · k_2,ω + 2 · k_3,ω + k_4,ω ) = 1.83 s^-1 + 0.02 s6 · ( 0.114 s^-2 + 2 · 0.293 s^-2 + 2 · 0.294 s^-2 + 0.474 s^-2 ) = 1.84 s^-1 < /p >

Kinetic energy - E_k ( i ) = m2 · l² · ⃗ω_i²

Potential energy- E_p ( i ) = m · l · g · ( 1 − cos ( ⃗θ_i ) )

...

Plot results

Rotation - θ(t), [°]

Energy - E(t), [J]

━━ Potential - E_p, ━━ Kinetic - E_k, ━━ Total - E_tot

Animate results

Simulation time - t = 0.14 s

Simulation time - t = 0.30 s

Simulation time - t = 0.46 s

Simulation time - t = 0.62 s

Simulation time - t = 0.78 s

Simulation time - t = 0.94 s

Simulation time - t = 1.10 s

Simulation time - t = 1.26 s

Simulation time - t = 1.42 s

Simulation time - t = 1.58 s

Simulation time - t = 1.74 s

Simulation time - t = 1.90 s

Simulation time - t = 2.06 s

Simulation time - t = 2.22 s

Simulation time - t = 2.38 s

Simulation time - t = 2.54 s

Simulation time - t = 2.70 s

Simulation time - t = 2.86 s

Simulation time - t = 3.02 s

Simulation time - t = 3.18 s

Simulation time - t = 3.34 s

Simulation time - t = 3.50 s

Simulation time - t = 3.66 s

Simulation time - t = 3.82 s

Simulation time - t = 3.98 s

Simulation time - t = 4.14 s

Simulation time - t = 4.30 s

Simulation time - t = 4.46 s

Simulation time - t = 4.62 s

Simulation time - t = 4.78 s

Simulation time - t = 4.94 s

Simulation time - t = 5.10 s

Simulation time - t = 5.26 s

Simulation time - t = 5.42 s

Simulation time - t = 5.58 s

Simulation time - t = 5.74 s

Simulation time - t = 5.90 s

Simulation time - t = 6.06 s

Simulation time - t = 6.22 s

Simulation time - t = 6.38 s

Simulation time - t = 6.54 s

Simulation time - t = 6.70 s

Simulation time - t = 6.86 s

Simulation time - t = 7.02 s

Simulation time - t = 7.18 s

Simulation time - t = 7.34 s

Simulation time - t = 7.50 s

Simulation time - t = 7.66 s

Simulation time - t = 7.82 s

Simulation time - t = 7.98 s

Simulation time - t = 8.14 s

Simulation time - t = 8.30 s

Simulation time - t = 8.46 s

Simulation time - t = 8.62 s

Simulation time - t = 8.78 s

Simulation time - t = 8.94 s

Simulation time - t = 9.10 s

Simulation time - t = 9.26 s

Simulation time - t = 9.42 s

Simulation time - t = 9.58 s

Simulation time - t = 9.74 s

Simulation time - t = 9.90 s

Simulation time - t = 10.06 s

Simulation time - t = 10.22 s

Simulation time - t = 10.38 s

Simulation time - t = 10.54 s

Simulation time - t = 10.70 s

Simulation time - t = 10.86 s

Simulation time - t = 11.02 s

Simulation time - t = 11.18 s

Simulation time - t = 11.34 s

Simulation time - t = 11.50 s

Simulation time - t = 11.66 s

Simulation time - t = 11.82 s

Simulation time - t = 11.98 s

Simulation time - t = 12.14 s

Simulation time - t = 12.30 s

Simulation time - t = 12.46 s

Simulation time - t = 12.62 s

Simulation time - t = 12.78 s

Simulation time - t = 12.94 s

Simulation time - t = 13.10 s

Simulation time - t = 13.26 s

Simulation time - t = 13.42 s

Simulation time - t = 13.58 s

Simulation time - t = 13.74 s

Simulation time - t = 13.90 s

Simulation time - t = 14.06 s

Simulation time - t = 14.22 s

Simulation time - t = 14.38 s

Simulation time - t = 14.54 s

Simulation time - t = 14.70 s

Simulation time - t = 14.86 s

Simulation time - t = 15.02 s

Simulation time - t = 15.18 s

Simulation time - t = 15.34 s

Simulation time - t = 15.50 s

Simulation time - t = 15.66 s

Simulation time - t = 15.82 s

Simulation time - t = 15.98 s

Simulation time - t = 16.14 s

Simulation time - t = 16.30 s

Simulation time - t = 16.46 s

Simulation time - t = 16.62 s

Simulation time - t = 16.78 s

Simulation time - t = 16.94 s

Simulation time - t = 17.10 s

Simulation time - t = 17.26 s

Simulation time - t = 17.42 s

Simulation time - t = 17.58 s

Simulation time - t = 17.74 s

Simulation time - t = 17.90 s

Simulation time - t = 18.06 s

Simulation time - t = 18.22 s

Simulation time - t = 18.38 s

Simulation time - t = 18.54 s

Simulation time - t = 18.70 s

Simulation time - t = 18.86 s

Simulation time - t = 19.02 s

Simulation time - t = 19.18 s

Simulation time - t = 19.34 s

Simulation time - t = 19.50 s

Simulation time - t = 19.66 s

Simulation time - t = 19.82 s

Simulation time - t = 19.98 s

Elastic Damped Pendulum 🎬¶

Two-DOF elastic pendulum: the rigid rod is replaced by an axial spring of stiffness $k_s$, giving coupled radial and angular equations of motion with independent damping factors. RK4 integration drives the SVG animation.

Code:

#rad
'<table><tr><td>
'<h4>Input parameters</h4>
'Gravitational acceleration -'g = 9.80665m/s^2
'Pendulum unstretched length -'l = 1m
'Pendulum mass -'m = 1kg
'Strut stiffness -'k_s = 40*(N/.m)
'Initial angle -'θ_0 = -179.9°|rad
'Initial angular velocity -'ω_0 = 0rad/s
'Initial length -'r_0 = 1.1*l
'Initial radial velocity -'v_0 = 0m/s
'Damping factor (normalized, radial) -'β_r = 0.05
'Damping factor (normalized, angular) -'β_θ = 0.01
'</td><td>
#hide
PlotWidth = 220','PlotHeight = 140','PlotSVG = 1
w = 1.2*l', 'h_ = 2/3*w
W = 150', 'H = 150
k = W/w
L = r_0*k
R = L/12
x = L*sin(θ_0)
y = L*cos(θ_0)
δ = 16
δ_x = 0.25*δ
δ_y = -2*δ
y_g = y - δ_y + R
#def svg$ = '<svg viewbox="'-100 - W/2' '-H' 'W + 200' '2*H  + 100'" xmlns="http://www.w3.org/2000/svg" version="1.1" style="font-family: Georgia Pro; font-size:14px; width:'W + 100'pt; height:'H + 100'pt;margin-left:-10pt;">
#def thin_style$ = style="stroke:goldenrod; stroke-width:1; fill:none;"
#def thick_style$ = style="stroke:steelblue; stroke-width:1.5; fill:none;"
#def ball_style$ = style="stroke:steelblue; stroke-width:0.5;  fill:url(#ball);"
#def bal0_style$ = style="stroke:steelblue; stroke-width:0.25;  fill:url(#ball); fill-opacity:0.5;"
#def solid_style$ = style="stroke:black; stroke-width:1; fill:#eee;"
#def load_style$ = style="stroke:steelblue; stroke-width:1; fill:steelblue;"
#include svg_drawing.cpd
#show
#val
'<br/><br/>
svg$
'<defs><radialGradient id="ball" cx="35%" cy="35%"><stop offset="0%" stop-color="lightcyan"/><stop offset="100%" stop-color="lightblue"/></radialGradient><pattern id="concrete" x="4" y="4" width="8" height="8" patternUnits="userSpaceOnUse"><rect x="0" y="0" width="8" height="8" fill="#eae9e8" /><circle cx="1" cy="1" r="1.2" fill="#ccb" /><circle cx="5" cy="2" r="1.6" fill="#dadad0" /><circle cx="4" cy="6" r="0.8" fill="#aa9" /><circle cx="3" cy="4" r="0.4" fill="#884" /><circle cx="7" cy="5" r="1.2" fill="#cacaba" /><circle cx="5" cy="3" r="0.9" fill="#fffded" /></pattern></defs>
'<polygon points="0,0 '-δ','-δ' 'δ','-δ'" solid_style$/>
line$(0; 0; 0; 0.4*L; thin_style$)
'<path d="M 0 '0.2*L' A '0.2*L' '0.2*L' 0 0 'θ_0 < 0°' '0.2*x' '0.2*y'" thin_style$/>
texth$(-1.5*δ*sign(θ_0);0.25*L;'θ_0*180/π'°)
line$(0; 0; x; y; thick_style$)
text$(x/2-δ_x*y/L; y/2+δ_x*x/L;450°+θ_0;'r_0' ('l') m)
circle$(0; 0; R/2; solid_style$)
circle$(0; 0; R/8; solid_style$)
circle$(l*k*sin(θ_0); l*k*cos(θ_0); R; bal0_style$)
circle$(x; y; R; ball_style$)
'<polygon points="'x','y_g + δ_y' 'x','y_g' 'x - δ_x','y_g + δ_y/3' 'x + δ_x','y_g + δ_y/3' 'x','y_g'" load_style$/>
texth$(x+3*δ_x; y_g-δ;g)
'</svg>
#equ
'</td></tr></table>
'Critical damping (radial) -'b_cr_r = 2*m*sqrt(k_s/m)
'Critical damping (angular) -'b_cr_θ = 2*m*l*sqrt(g/l)
'Radial damping -'b_r = β_r*b_cr_r
'Angular damping -'b_θ = β_θ*b_cr_θ
'Maximum simulation time -'t_max = 20s
'<h4>Solution by Runge-Kutta RK4 method (explicit)</h4>
'Differential equations
'<div style="border-left:solid 1pt black; margin-left:24pt; padding-left:3pt;">
#noc
m*r″ - m*r*θ′^2 + k_s*(r - l) + b_r*r′ - m*g*cos(θ) = 0
m*r*θ″ + 2*m*r′*θ′ + b_θ*θ′ + m*g*sin(θ) = 0
#equ
'</div>
'The II order system of 2 equations is first reduced to the following I order system of 4 equations
'<div style="border-left:solid 1pt black; margin-left:24pt; padding-left:3pt;">
#noc
r′ = v
v′ = r*ω^2 - k/m*(r - l) - b_r/m*v + g*cos(θ)
θ′ = ω
ω′ = -((2*m*v*ω + b_θ*ω + m*g*sin(θ))/(m*r))
#equ
'</div>
'The solution is performed iteratively
'Step size -'h = 0.01s
'Number of steps -'n = t_max/h
#hide
N = min(n; 120)
k_N = floor(n/N)
N = n\k_N
#show
'Allocate vectors
r = vector(n)
v = vector(n)
θ = vector(n)
ω = vector(n)
'Set initial conditions
r.1 = r_0','θ.1 = θ_0/1rad
v.1 = v_0','ω.1 = ω_0/1rad
'Define the system of ODEs
f_r(r; v; θ; ω) = v
f_v(r; v; θ; ω) = r*ω^2 + g*cos(θ) - k_s/m*(r - l) - b_r/m*v
f_θ(r; v; θ; ω) = ω
f_ω(r; v; θ; ω) = -((2*m*v*ω + b_θ*ω + m*g*sin(θ))/(m*r))
'Perform Runge-Kutta 4 steps
#for i = 1 : n - 1
    r_i = r.i', 'v_i = v.i
    θ_i = θ.i', 'ω_i = ω.i
    k_1r = h*v_i', 'k_1v = h*f_v(r_i; v_i; θ_i; ω_i)
    k_1θ = h*ω_i', 'k_1ω = h*f_ω(r_i; v_i; θ_i; ω_i)
    r_i = r.i + k_1r/2', 'v_i = v.i + k_1v/2
    θ_i = θ.i + k_1θ/2', 'ω_i = ω.i + k_1ω/2
    k_2r = h*v_i', 'k_2v = h*f_v(r_i; v_i; θ_i; ω_i)
    k_2θ = h*ω_i', 'k_2ω = h*f_ω(r_i; v_i; θ_i; ω_i)
    r_i = r.i + k_2r/2', 'v_i = v.i + k_2v/2
    θ_i = θ.i + k_2θ/2', 'ω_i = ω.i + k_2ω/2
    k_3r = h*v_i', 'k_3v = h*f_v(r_i; v_i; θ_i; ω_i)
    k_3θ = h*ω_i', 'k_3ω = h*f_ω(r_i; v_i; θ_i; ω_i)
    r_i = r.i + k_3r', 'v_i = v.i + k_3v
    θ_i = θ.i + k_3θ', 'ω_i = ω.i + k_3ω
    k_4r = h*v_i', 'k_4v = h*f_v(r_i; v_i; θ_i; ω_i)
    k_4θ = h*ω_i', 'k_4ω = h*f_ω(r_i; v_i; θ_i; ω_i)
    'Update values for the next step
    r.(i + 1) = r.i + (k_1r + 2*k_2r + 2*k_3r + k_4r)/6
    v.(i + 1) = v.i + (k_1v + 2*k_2v + 2*k_3v + k_4v)/6
    θ.(i + 1) = θ.i + (k_1θ + 2*k_2θ + 2*k_3θ + k_4θ)/6
    ω.(i + 1) = ω.i + (k_1ω + 2*k_2ω + 2*k_3ω + k_4ω)/6
    #hide
#loop
t(i) = i*h - h
r2d = 180/π
#show
'Kinetic energy -'E_k(i) = m/2*(v.i^2 + (r.i*ω.i)^2)|J
'Potential energy -'E_p(i) = k_s/2*(r.i - l)^2 - m*g*(r.i*cos(θ.i) - l)|J
'...
'<h4>Plot results</h4><table><tr><td>Angle - <var>θ</var>(<var>t</var>), [°]<br/>
$Plot{t(i)|θ.i*r2d @ i = 1 : n}
'</td><td>Length - <var>r</var>(<var>t</var>), [m]<br/>
$Plot{t(i)|r.i @ i = 1 : n}
'</td></tr><tr><td>Trajectory - <var>x</var>(<var>t</var>), <var>y</var>(<var>t</var>), [m]<br/>
$Plot{r.i*sin(θ.i)|-r.i*cos(θ.i) @ i = 1 : n}
'</td><td>Energy - <var>E</var>(<var>t</var>), [J]<br/>
$Plot{t(i)|E_p(i) + E_k(i) & 0s|0J & 0s|0J & 0s|0J & t(i)|E_p(i) & t(i)|E_k(i) @ i = 1 : n}
'<br/><b style="color:MediumVioletRed">━━</b> Potential - <var>E</var><sub>p</sub>,
'<b style="color:MediumSpringGreen">━━</b> Kinetic - <var>E</var><sub>k</sub>,
'<b style="color:Tomato">━━</b> Total - <var>E</var><sub>tot</sub>
'</td></tr></table>
'<style>.Series3, .Series4{stroke:none!Important; fill:none!Important;} .fr{display:none;} .fr svg{align:center;} .trace{stroke:red; stroke-width:1; fill:red;}</style>
'<h4>Animate results</h4>
#format 0.00
#hide
k = W/w
x = r*sin(θ)
y = r*cos(θ)
n_tr = 40
#for i_N = 1 : N
    #hide
    i = i_N*k_N
    xi = x.i*k
    yi = y.i*k
    t = t(i)
    i_0 = max((i_N - n_tr)*k_N; 2)
    n_0 = i - i_0
    Δr = abs(r_0 - l) + 0.2*l
    #show
    #val
    '<div class="fr" id="p4-fr-'i_N'">
    #equ
    #nosub
    'Simulation time -'t'<br/>
    '<table><tr><td>Angle - <var>θ</var>(<var>t</var>), [°]<br/>
    $Plot{t(i)|θ.i*r2d & t(j)|θ.j*r2d & t_max|θ_0/1rad*r2d & t_max|-θ_0/1rad*r2d @ j = 1 : i}
    '<br/>Length - <var>r</var>(<var>t</var>), [m]<br/>
    $Plot{t(i)|r.i & t(j)|r.j & t_max|l + Δr & t_max|l - Δr @ j = 1 : i}
    '</td><td>
    #val
    svg$
    #for j = 1 : n_0/2
        #hide
        j2 = i_0 + 2*j
        j1 = j2 - 2
        o = round(j/n_0*100)/100
        #show
        '<line x1="'x.j1*k'" y1="'y.j1*k'" x2="'x.j2*k'" y2="'y.j2*k'" class="trace" opacity="'o'"/>
    #loop
    '<polygon points="0,0 '-δ','-δ' 'δ','-δ'" solid_style$/>
    line$(0; 0; xi; yi; thick_style$)
    circle$(0; 0; R/2; solid_style$)
    circle$(0; 0; R/8; solid_style$)
    circle$(xi; yi; R; ball_style$)
    '</svg></td></tr></table></div>
#loop
'<script>(function(){$("#p4-fr-1").show();var i=1;var fr=$("#p4-fr-1");setInterval(function(){fr.hide();if(++i>'N')i=1;fr=$("#p4-fr-"+i);fr.show();},'1000*h*k_N');})();</script>
#equ

Rendered Output:

Input parameters

Gravitational acceleration - g = 9.81 ms²

Pendulum unstretched length - l = 1 m

Pendulum mass - m = 1 kg

Strut stiffness - k_s = 40 Nm

Initial angle - θ₀ = -179.9° = -3.14 rad

Initial angular velocity - ω₀ = 0 rads

Initial length - r₀ = 1.1 · l = 1.1 · 1 m = 1.1 m

Initial radial velocity - v₀ = 0 ms

Damping factor (normalized, radial) - β_r = 0.05

Damping factor (normalized, angular) - β_θ = 0.01

Critical damping (radial) - b_{cr_r} = 2 · m · k_sm = 2 · 1 kg · 40 N ∕ m1 kg = 12.65 kg ∕ s

Critical damping (angular) - b_{cr_θ} = 2 · m · l · gl = 2 · 1 kg · 1 m · 9.81 m ∕ s²1 m = 6.26 kg · m ∕ s

Radial damping - b_r = β_r · b_{cr_r} = 0.05 · 12.65 kg ∕ s = 0.632 kg ∕ s

Angular damping - b_θ = β_θ · b_{cr_θ} = 0.01 · 6.26 kg · m ∕ s = 0.0626 kg · m ∕ s

Maximum simulation time - t_max = 20 s

Solution by Runge-Kutta RK4 method (explicit)

Differential equations

m · r″ − m · r · θ′² + k_s · ( r − l ) + b_r · r′ − m · g · cos ( θ ) = 0

m · r · θ″ + 2 · m · r′ · θ′ + b_θ · θ′ + m · g · sin ( θ ) = 0

The II order system of 2 equations is first reduced to the following I order system of 4 equations

r′ = v

v′ = r · ω² − km · ( r − l ) − b_rm · v + g · cos ( θ )

θ′ = ω

ω′ = - 2 · m · v · ω + b_θ · ω + m · g · sin ( θ ) m · r

The solution is performed iteratively

Step size - h = 0.01 s

Number of steps - n = t_maxh = 20 s0.01 s = 2000

Allocate vectors

⃗r = vector ( n ) = vector ( 2000 ) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... 0]

⃗v = vector ( n ) = vector ( 2000 ) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... 0]

⃗θ = vector ( n ) = vector ( 2000 ) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... 0]

⃗ω = vector ( n ) = vector ( 2000 ) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... 0]

Set initial conditions

⃗r₁ = r₀ = 1.1 m , ⃗θ₁ = θ₀1 rad = -3.14 rad1 rad = -3.14

⃗v₁ = v₀ = 0 m ∕ s , ⃗ω₁ = ω₀1 rad = 0 rad ∕ s1 rad = 0 s^-1

Define the system of ODEs

f_r ( r; v; θ; ω ) = v

f_v ( r; v; θ; ω ) = r · ω² + g · cos ( θ ) − k_sm · ( r − l ) − b_rm · v

f_θ ( r; v; θ; ω ) = ω

f_ω ( r; v; θ; ω ) = - 2 · m · v · ω + b_θ · ω + m · g · sin ( θ ) m · r

Perform Runge-Kutta 4 steps

r_i = ⃗r₁ = 1.1 m , v_i = ⃗v₁ = 0 m ∕ s

θ_i = ⃗θ₁ = -3.14 , ω_i = ⃗ω₁ = 0 s^-1

k_1r = h · v_i = 0.01 s · 0 m ∕ s = 0 m , k_1v = h · f_v ( r_i; v_i; θ_i; ω_i ) = 0.01 s · f_v ( 1.1 m; 0 m ∕ s; -3.14; 0 s^-1 ) = -0.138 m ∕ s

k_1θ = h · ω_i = 0.01 s · 0 s^-1 = 0 , k_1ω = h · f_ω ( r_i; v_i; θ_i; ω_i ) = 0.01 s · f_ω ( 1.1 m; 0 m ∕ s; -3.14; 0 s^-1 ) = 0.000156 s^-1

r_i = ⃗r₁ + k_1r2 = 1.1 m + 0 m2 = 1.1 m , v_i = ⃗v₁ + k_1v2 = 0 m ∕ s + -0.138 m ∕ s2 = -0.069 m ∕ s

θ_i = ⃗θ₁ + k_1θ2 = -3.14 + 02 = -3.14 , ω_i = ⃗ω₁ + k_1ω2 = 0 s^-1 + 0.000156 s^-12 = 7.78×10^-5 s^-1

k_2r = h · v_i = 0.01 s · ( -0.069 m ∕ s ) = -0.00069 m , k_2v = h · f_v ( r_i; v_i; θ_i; ω_i ) = 0.01 s · f_v ( 1.1 m; -0.069 m ∕ s; -3.14; 7.78×10^-5 s^-1 ) = -0.138 m ∕ s

k_2θ = h · ω_i = 0.01 s · 7.78×10^-5 s^-1 = 7.78×10^-7 , k_2ω = h · f_ω ( r_i; v_i; θ_i; ω_i ) = 0.01 s · f_ω ( 1.1 m; -0.069 m ∕ s; -3.14; 7.78×10^-5 s^-1 ) = 0.000156 s^-1

r_i = ⃗r₁ + k_2r2 = 1.1 m + -0.00069 m2 = 1.1 m , v_i = ⃗v₁ + k_2v2 = 0 m ∕ s + -0.138 m ∕ s2 = -0.0688 m ∕ s

θ_i = ⃗θ₁ + k_2θ2 = -3.14 + 7.78×10^-72 = -3.14 , ω_i = ⃗ω₁ + k_2ω2 = 0 s^-1 + 0.000156 s^-12 = 7.78×10^-5 s^-1

k_3r = h · v_i = 0.01 s · ( -0.0688 m ∕ s ) = -0.000688 m , k_3v = h · f_v ( r_i; v_i; θ_i; ω_i ) = 0.01 s · f_v ( 1.1 m; -0.0688 m ∕ s; -3.14; 7.78×10^-5 s^-1 ) = -0.137 m ∕ s

k_3θ = h · ω_i = 0.01 s · 7.78×10^-5 s^-1 = 7.78×10^-7 , k_3ω = h · f_ω ( r_i; v_i; θ_i; ω_i ) = 0.01 s · f_ω ( 1.1 m; -0.0688 m ∕ s; -3.14; 7.78×10^-5 s^-1 ) = 0.000156 s^-1

r_i = ⃗r₁ + k_3r = 1.1 m + ( -0.000688 m ) = 1.1 m , v_i = ⃗v₁ + k_3v = 0 m ∕ s + ( -0.137 m ∕ s ) = -0.137 m ∕ s

θ_i = ⃗θ₁ + k_3θ = -3.14 + 7.78×10^-7 = -3.14 , ω_i = ⃗ω₁ + k_3ω = 0 s^-1 + 0.000156 s^-1 = 0.000156 s^-1

k_4r = h · v_i = 0.01 s · ( -0.137 m ∕ s ) = -0.00137 m , k_4v = h · f_v ( r_i; v_i; θ_i; ω_i ) = 0.01 s · f_v ( 1.1 m; -0.137 m ∕ s; -3.14; 0.000156 s^-1 ) = -0.137 m ∕ s

k_4θ = h · ω_i = 0.01 s · 0.000156 s^-1 = 1.56×10^-6 , k_4ω = h · f_ω ( r_i; v_i; θ_i; ω_i ) = 0.01 s · f_ω ( 1.1 m; -0.137 m ∕ s; -3.14; 0.000156 s^-1 ) = 0.000156 s^-1

Update values for the next step

⃗r₂ = ⃗r₁ + k_1r + 2 · k_2r + 2 · k_3r + k_4r6 = 1.1 m + 0 m + 2 · ( -0.00069 m ) + 2 · ( -0.000688 m ) + ( -0.00137 m ) 6 = 1.1 m

⃗v₂ = ⃗v₁ + k_1v + 2 · k_2v + 2 · k_3v + k_4v6 = 0 m ∕ s + -0.138 m ∕ s + 2 · ( -0.138 m ∕ s ) + 2 · ( -0.137 m ∕ s ) + ( -0.137 m ∕ s ) 6 = -0.138 m ∕ s

⃗θ₂ = ⃗θ₁ + k_1θ + 2 · k_2θ + 2 · k_3θ + k_4θ6 = -3.14 + 0 + 2 · 7.78×10^-7 + 2 · 7.78×10^-7 + 1.56×10^-66 = -3.14

⃗ω₂ = ⃗ω₁ + k_1ω + 2 · k_2ω + 2 · k_3ω + k_4ω6 = 0 s^-1 + 0.000156 s^-1 + 2 · 0.000156 s^-1 + 2 · 0.000156 s^-1 + 0.000156 s^-16 = 0.000156 s^-1

Kinetic energy - E_k ( i ) = m2 · ( ⃗v_i² + ( ⃗r_i · ⃗ω_i ) ² )

Potential energy - E_p ( i ) = k_s2 · ( ⃗r_i − l ) ² − m · g · ( ⃗r_i · cos ( ⃗θ_i ) − l )

...

Plot results

Angle - `θ`(`t`), [°]	Length - `r`(`t`), [m]
Trajectory - `x`(`t`), `y`(`t`), [m]	Energy - `E`(`t`), [J] ━━ Potential - `E`_p, ━━ Kinetic - `E`_k, ━━ Total - `E`_tot

Animate results

Simulation time - t = 0.15 s