- 624

- 142

No external dependencies.

Arbitrary order simulation (accuracy limited by float precision).

No finite difference errors.

Can be extended to arbitrary precision with gmpy2 and two more lines.

Enjoy!

Parameters:

- order
- step size
- number of steps
- initial conditions (x, y, z)
- parameters (sigma, rho, beta_numerator, beta_denominator)

Code:

`./tsm-lorenz.py 10 .01 10001 -15.8 -17.48 35.64 10 28 8 3`

Python:

```
#!/usr/bin/env python3
from sys import argv
n = int(argv[1])
h = float(argv[2])
steps = int(argv[3])
x = float(argv[4])
y = float(argv[5])
z = float(argv[6])
s = float(argv[7])
r = float(argv[8])
b = float(argv[9]) / float(argv[10])
cx = [0.0 for _ in range(n + 1)]
cy = [0.0 for _ in range(n + 1)]
cz = [0.0 for _ in range(n + 1)]
for t in range(1, steps):
print("{:.9e} {:.9e} {:.9e} {:.5e}".format(x, y, z, t * h))
cx[0] = x
cy[0] = y
cz[0] = z
for k in range(n):
xz = xy = 0.0
cx[k + 1] = s * (cy[k] - cx[k]) / (k + 1)
for j in range(k + 1):
xz += cx[k - j] * cz[j]
cy[k + 1] = (r * cx[k] - xz - cy[k]) / (k + 1)
for j in range(k + 1):
xy += cx[k - j] * cy[j]
cz[k + 1] = (xy - b * cz[k]) / (k + 1)
x = cx[n]
y = cy[n]
z = cz[n]
for i in range(n - 1, -1, -1):
x = x * h + cx[i]
y = y * h + cy[i]
z = z * h + cz[i]
```