The Wimol-Banlue attractor via the Taylor Series Method

[m4r35n357]

This attractor is unusual because it uses both the tanh() and abs() functions. A picture can be found here (penultimate image). Here is some dependency-free Python (abridged from the GitHub code, but not flattened!) to generate the data to arbitrary order:

#!/usr/bin/env python3
from sys import argv
from math import tanh

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

def t_abs(u, k):
    if k == 0:
        return abs(u[0])
    elif k == 1:
        return u[1] if u[0] > 0.0 else (-u[1] if u[0] < 0.0 else 0.0)
    return 0.0

def t_prod(u, v, k):
    return sum(u[j] * v[k - j] for j in range(k + 1))

def t_tanh_sech2(t, s2, u, k):
    if k == 0:
        return tanh(u[0]), 1.0 - tanh(u[0])**2
    tk = s2[0] * u[k] + sum(j * u[j] * s2[k - j] for j in range(1, k)) / k
    return tk, - (2.0 * (t[0] * tk + sum(j * t[j] * t[k - j] for j in range(1, k)) / k))

def main():
    model, order, δt, n_steps = argv[1], int(argv[3]), float(argv[4]), int(argv[5])  # integrator controls
    x0, y0, z0 = float(argv[6]), float(argv[7]), float(argv[8])  # initial values
    x, y, z = t_jet(order + 1), t_jet(order + 1), t_jet(order + 1)  # coordinate jets

    if model == "wimol-banlue":
            a_ = t_jet(order, float(argv[9]))  # constant jet
            tx, sx = t_jet(order), t_jet(order)  # temporary jets
            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  # Taylor Series Method
                for k in range(order):
                    tx[k], sx[k] = t_tanh_sech2(tx, sx, x, k)  # sech^2 stored but not used in ODE equations
                    x[k + 1] = (y[k] - x[k]) / (k + 1)
                    y[k + 1] = - t_prod(z, tx, k) / (k + 1)
                    z[k + 1] = (- a_[k] + t_prod(x, y, k) + t_abs(y, k)) / (k + 1)
                x0, y0, z0 = x[order], y[order], z[order]  # Horner's method
                for i in range(order - 1, -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}')

main()

To run it

./wb.py wimol-banlue 16 10 0.1 10001 1 0 0 2.0

Enjoy!
[EDIT] fixed bug in t_abs() !

 

[jedishrfu]

Very interesting attractors.

Can you provide a sample arg list to try?

I tried it on my iPad using Pythonista hoping it might run but the args stopped me.

 

[jedishrfu]

Also it doesn’t plot anything I guess that part. is left to the student.

 

[m4r35n357]

Can you provide a sample arg list to try?

See OP ;)
Also, the posters on the website in the OP have parameters on them.

 

[m4r35n357]

Also it doesn’t plot anything I guess that part. is left to the student.

You expect much from 45 lines of Python ;)
I have plotting programs on GitHub (there is a branch called float if you don't want arbitrary precision), but will you need matplotlib for graphs and Pi3d for real time 3d trajectories (tentative virtual env instructions in the README).
Otherwise, the data should be easy enough to parse into something else for display (SciLab?).

 

[jedishrfu]

Okay thanks. I use python for work. Mostly the Anaconda distro although I also have Pythonista on the iPad. Both have numpy and matplotlib modules so it should be possible to generate a plot which would make the code even more interesting visually.

 

[m4r35n357]

Okay thanks. I use python for work. Mostly the Anaconda distro although I also have Pythonista on the iPad. Both have numpy and matplotlib modules so it should be possible to generate a plot which would make the code even more interesting visually.

OK, so if you have matplotlib you should at least be able to plot the graphs in progressive animated fashion using the plotAnimated.py script. See the README for an example invocation.
pi3d is much more impressive visually though.