拡散項を含む1次元複素TDGL方程式の描画pyhthon script

描画にmatplotlib, パラメータ指定等のGUIにtkinterを使っています。

#!/usr/bin/python3
'''
Calculation and Animation for One-dimensional complex TDGL equation
(solved by  simple Euler method.)

'''
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import animation
import sys
import math


size = 400
amp = np.zeros(size)    # used in the computing process as amplitudes
diff = np.array(amp, dtype=np.complex128)   # used as the Laplacian term

x = np.random.rand(size)
x = np.array(x, dtype=np.complex128)
y = np.array(x, dtype=np.complex128)
site_list = range(size)
# parameters of the differential equation
c0 = 1.0
c1 = -1.0
c2 = 0.6
c0i = 1.0 + c0*1.0j
Di = 1.0  + c1*1.0j
ci = 1.0 + c2*1.0j
boundary = "f"

def updateField(x_in, x_out, c0i =0.0 + 0.0j,  ci=1 + 0.6j, Di= 1 -1.0j, boundary="f"):
    DT = 0.002 # time step
    amp = np.absolute(x_in)
    amp = np.power(amp,2)
    diff[1:-1] = x_in[0:-2,] + x_in[2:]  - 2.0* x_in[1:-1]
    x_out[1:-1] = x_in[1:-1] + DT*(c0i - ci* amp[1:-1])*x_in[1:-1] + DT* Di * diff[1:-1]
    # periodic boundary condition
    if boundary=="periodic":
        x_out[0] = x_out[-2]
        x_out[-1] = x_out[1]
    elif boundary == "fixed":
        x_out[0] =0.2
        x_out[-1] = 0.2
    else:
        x_out[0] = 0.2
        x_out[-1] = x_out[-2]

fig = plt.figure(figsize=(7,4))
ax = fig.add_subplot(1,1,1)
ax.set_ylim(-2.0, 2.0)
plt_real_x, = ax.plot(site_list, np.real(x))

# for interpolation, see  https://matplotlib.org/gallery/images_contours_and_fields/interpolation_methods.html
# for cmap,

def animate(i):
    global x, y, c0i, ci, Di, size, boundary

    # integration
    for count in range(20):
        updateField(x, y, c0i, ci, Di, boundary)
        updateField(y, x, c0i,  ci, Di, boundary)

    # update data for plotting
    # crgb =  (np.angle(x) + math.pi)/2.0/math.pi
    plt_real_x.set_data(site_list, np.real(x))
    return [plt_real_x]


if __name__ == '__main__':
    args = sys.argv
    if len(args) < 6:
        print("Usage: ctdgl_pde2d.py c0 c1 c2 boundary initial_condition")
        print(" typical values 1. 0.0 0.6 free 0")
        print(' boundary is "free" or "fixed" of "pediodic" ')
        print(' initial_condition 0 or 1, which means zeros or random')
        sys.exit()
    c0 = float(args[1]) # -0.5
    c1 = float(args[2]) # 0.0
    c2 = float(args[3]) # 0.6
    boundary = args[4]
    initial_condition = int(args[5])
    if not boundary in ["free", "fixed", "periodic"]:
        print('boundary must be one of "free", "fixed", "periodic"')
        sys.exit()

    c0i = 1.0 + c0*1.0j
    Di = 1.0  + c1*1.0j
    ci = 1.0 + c2*1.0j
    plt.title("$c_0=" + str(c0) + ", c_1=" + str(c1) + ", c_2=" + str(c2) + "$, boudary=" + boundary)
    # initial configuration
    if initial_condition == 1:
        for n in range(size):
            x[n] = 0.01*(np.random.rand() -0.5 + (np.random.random() - 0.5)* 1.0j)
    else:
        for n in range(size):
            x[n] = 0.0 + 0.0j

    anim = animation.FuncAnimation(fig, animate,
                               frames=1000, interval=10)
    # anim.save('ctdgl2d.c1_1.0_c2_0.6.mp4', fps=15, extra_args=['-vcodec', 'libx264'])
    # anim.save('oscillator_xyz_t.gif', writer='imagemagick',  fps=30)
    plt.show()