A Challenge to RK4

Hacker News

Share card

A Challenge to RK4

Why limit yourself to fourth-order integration? Why put up with roundoff errors in your finite differences? Here is a Lorenz ODE solver in 20 lines of dependency-free Python to enlighten you! Find out more about the Taylor Series Method at https://github.com/m4r35n357/ODE-Playground, read the docs and follow the references, marvel at how long this technique has lain undiscovered by the masses. #!/usr/bin/env python3 from sys import argv # Example: ./tsm.py lorenz 10 .01 3000 -15.8 -17.48 35.64 10 28 8 3 model, order, δt, n_steps = argv[1], int(argv[2]), float(argv[3]), int(argv[4]) # integrator controls x0, y0, z0 = float(argv[5]), float(argv[6]), float(argv[7]) # initial values x, y, z = [0.0] * (order + 1), [0.0] * (order + 1), [0.0] * (order + 1) # Taylor coefficient jets σ, ρ, β = float(argv[8]), float(argv[9]), float(argv[10]) / float(argv[11]) # Lorenz parameters print(f'{x0:.9e} {y0:.9e} {z0:.9e} {0.0:.5e}') for step in range(1, n_steps + 1): x[0], y[0], z[0] = x0, y0, z0 # build the coefficient jets for k in range(order): x[k + 1] = σ * (y[k] - x[k]) / (k + 1) y[k + 1] = (ρ * x[k] - sum(x[j] * z[k - j] for j in range(k + 1)) - y[k]) / (k + 1) z[k + 1] = (sum(x[j] * y[k - j] for j in range(k + 1)) - β * z[k]) / (k + 1) x0, y0, z0 = x[-1], y[-1], z[-1] # perform an integration step using Horner's method for i in range(len(x) - 2, -1, -1): x0 = x0 * δt + x[i] y0 = y0 * δt + y[i] z0 = z0 * δt + z[i] print(f'{x0:.9e} {y0:.9e} {z0:.9e} {step * δt:.5e}')