Python+NumbaでWikipediaのList of chaotic mapsを１つずつ描く（その4） Ikeda Map

今日は

https://en.wikipedia.org/wiki/List_of_chaotic_maps

の中の、Ikeda Mapです。

https://en.wikipedia.org/wiki/Ikeda_map

t[i]=0.4-6./(1+x[i]**2+y[i]**2)
x[i+1]=1+u*(x[i]*cos(t[i])-y[i]*sin(t[i]))
y[i+1]=u*(x[i]*sin(t[i])+y[i]*cos(t[i]))

の形式を使って計算する。Wikipediaに乗っているu=0.918で初期値を振ってみた。

こっちは狭い領域で有名な形。

これは広い領域でWikipediaのトップの絵に相当するもの。

リストはこちら。

import numpy as np
import matplotlib.pyplot as plt
from numba import jit

up=0.918

@jit
def ikeda_calc(x,y,N):

    for i in range(N-1):
        t=0.4-6./(1+x[i]**2+y[i]**2)
        x[i+1]=1+up*(x[i]*np.cos(t)-y[i]*np.sin(t))
        y[i+1]=up*(x[i]*np.sin(t)+y[i]*np.cos(t))

N=500
xmax=10
xmin=-10
ymax=10
ymin=-10
Nx=100
Ny=100

fig=plt.figure(figsize=(10,10))
plt.xlim(-10,10)
plt.ylim(-10,10)

for i in range(Nx):
    for j in range(Ny):
        x=np.zeros(N)
        y=np.zeros(N)
        x[0],y[0]=xmin+(xmax-xmin)*i/Nx, ymin+(ymax-ymin)*j/Ny
        ikeda_calc(x,y,N)
        plt.scatter(x,y,s=1,c='black',alpha=0.01)

plt.show()