Damped & Driven Pendulums (in _pure_ Python)

  • #1
654
148

Summary:

Minimal example code for experimentation

Main Question or Discussion Point

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!

Simulator:
#!/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
 
  • Like
Likes Filip Larsen

Answers and Replies

  • #2
11,809
5,434
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##
 
  • #3
654
148
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:
Python:
#!/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]
 
Last edited:

Related Threads on Damped & Driven Pendulums (in _pure_ Python)

  • Last Post
Replies
0
Views
4K
Replies
9
Views
3K
Replies
17
Views
1K
Replies
6
Views
4K
  • Last Post
Replies
10
Views
1K
  • Last Post
Replies
4
Views
420
  • Last Post
Replies
22
Views
5K
Replies
7
Views
1K
  • Last Post
Replies
4
Views
6K
  • Last Post
Replies
1
Views
2K
Top