- #1
m4r35n357
- 654
- 148
- 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!
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