Damped & Driven Pendulums (in _pure_ Python)

[m4r35n357]

TL;DR Summary

Minimal example code for experimentation

This is another application of using Taylor recurrences (open access) to solve ODEs to arbitrarily high order (e.g. 10th order in the example invocation). It illustrates use of trigonometric recurrences, rather than the product recurrences in my earlier Lorenz ODE posts.

Enjoy!

#!/usr/bin/env python3
#  Example: ./forced.py 10 .05 10001 1.0 1.0 0.1 4.9 1.1
from sys import argv
from math import pi, sqrt, sin, cos

def jet(n, value=0.0):
    j = [0.0] * n
    j[0] = value
    return j

def Σ(j, n, h):
    result = j[n]
    for i in range(n - 1, -1, -1):
        result = result * h + j[i]
    return result

def _dd(v, u, i):
    return sum(j * u[j] * v[i - j] for j in range(1, i)) / i

def t_sin_cos(s, c, u, i):
    return (sin(u[0]), cos(u[0])) if i == 0 else (c[0] * u[i] + _dd(c, u, i), - s[0] * u[i] - _dd(s, u, i))

order, δt, n_steps = int(argv[1]), float(argv[2]), int(argv[3])  # integrator controls
g, m, length = 9.80665, float(argv[4]), float(argv[5])  # physical parameters
ζ, a, ω = float(argv[6]), float(argv[7]), 2.0 * pi * sqrt(length / g) * float(argv[8])  # damping/forcing parameters

θ0, θdot0 = 0.0, 0.0  # initial values
θ, θdot, sinθ, cosθ = jet(order + 1), jet(order + 1), jet(order), jet(order)  # jets
for step in range(1, n_steps):
    print("{:.9e} {:.9e} {:.9e} {:.5e}".format(θ0, θdot0, float(0.0), step * δt))
    θ[0], θdot[0] = θ0, θdot0
    for k in range(order):  # build up jets using recurrences and the derivative rule
        sinθ[k], cosθ[k] = t_sin_cos(sinθ, cosθ, θ, k)
        θ[k + 1] = θdot[k] / (k + 1)
        θdot[k + 1] = (a * cos(ω * step * δt) - ζ * length * θdot[k] - m * g * sinθ[k]) / (m * length) / (k + 1)
    θ0, θdot0 = Σ(θ, order, δt), Σ(θdot, order, δt)  # Horner's method

 

[jedishrfu]

Thanks for sharing!

It’s always great to code snippets. I didn’t realize that in python you could use Greek letters for method and variable names like $\Sigma$

 

[m4r35n357]

Thanks for sharing!

Cheers. Here is a slightly tidied and reduced version, single-use functions are eliminated, and now uses f-strings for output formatting:

#!/usr/bin/env python3
#  Example: ./forced.py 10 .05 4001 1.0 1.0 0.1 4.9 1.1 | ./plotAnimated.py 1 -50 50
from sys import argv
from math import pi, sqrt, sin, cos

def dd(v, u, i):
    return sum(j * u[j] * v[i - j] for j in range(1, i)) / i

order, δt, n_steps = int(argv[1]), float(argv[2]), int(argv[3])  # integrator controls
g, m, l = 9.80665, float(argv[4]), float(argv[5])  # physical parameters
ζ, a, ω = float(argv[6]), float(argv[7]), 2.0 * pi * sqrt(l / g) * float(argv[8])  # damping/forcing

θ0, θdot0 = 0.0, 0.0  # initial values
θ, θdot, sinθ, cosθ = [0.0] * (order + 1), [0.0] * (order + 1), [0.0] * order, [0.0] * order  # jets
for step in range(1, n_steps):
    print(f'{θ0:.9e} {θdot0:.9e} {0.0:.9e} {step * δt:.5e}')

    # build up jets using recurrences and the derivative rule
    θ[0], θdot[0] = θ0, θdot0
    for k in range(order):
        sinθ[k], cosθ[k] = (sin(θ[0]), cos(θ[0])) if k == 0 else (cosθ[0] * θ[k] + dd(cosθ, θ, k),
                                                                - sinθ[0] * θ[k] - dd(sinθ, θ, k))
        θ[k + 1] = θdot[k] / (k + 1)
        θdot[k + 1] = (a * cos(ω * step * δt) - ζ * l * θdot[k] - m * g * sinθ[k]) / (m * l) / (k + 1)

    # evaluate series using Horner's method
    θ0, θdot0 = θ[order], θdot[order]
    for k in range(order - 1, -1, -1):
        θ0, θdot0 = θ0 * δt + θ[k], θdot0 * δt + θdot[k]