Simulation of the Time-Dependent Ginzburg-Landau Equation $$\frac{\text{d}u}{\text{d}t}= u - u^3 +D\nabla^2 u,$$ in 1 and 2 spatial dimensions. This is the simplest example of numerical integration through Finite Differences:

Written by Yair Mau. Check out my webpage for more tutorials: http://www.yairmau.com/

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import time as tm

In [2]:
# grid size
n = 128
# for 1d simulation write N=(n,)
N=(n,n)
# diffusion coefficient
D = 1.0
# spatial dimensions
L = 100.0
dx = L / n
x = np.arange(0,L,dx)
# time
t = 0.0
total_time = 3.0
# beware of the Von Neumann stability analysis
# https://en.wikipedia.org/wiki/Von_Neumann_stability_analysis
dt = 0.2 * 0.5 * dx**2 / D

In [3]:
# define functions
def periodic_lap_1d(u,dx=1.0):
return (+1*np.roll(u,+1)
+1*np.roll(u,-1)
-2*u) / dx**2
def periodic_lap_2d(u,dx=1.0):
return (+1*np.roll(u,+1,axis=0)
+1*np.roll(u,-1,axis=0)
+1*np.roll(u,+1,axis=1)
+1*np.roll(u,-1,axis=1)
-4*u) / dx**2
f = lambda u: u - u**3

In [4]:
# initialize and start plotting
plt.ion()
fig = plt.figure(1,figsize=(7,6))
plt.clf()
# random initial condition
u = 2*np.random.random(N)-1.0
if len(N) == 1:
lap = periodic_lap_1d
p, = ax.plot(x,u)
ax.axis([x[0],x[-1],-1.1,1.1])
if len(N) == 2:
lap = periodic_lap_2d
p = ax.imshow(u,cmap="RdGy", vmin=-1.0, vmax=1.0,extent=[0,L,0,L])
# create an axes on the right side of ax. The width of cax will be 5%
# of ax and the padding between cax and ax will be fixed at 0.15 inch.
divider = make_axes_locatable(ax)
cbar = fig.colorbar(p, cax=colorbar_ax, ticks=[-1,-0.5,0,0.5,1])
ax.set_title("time={:5.1f}".format(0.0))

Out[4]:
<matplotlib.text.Text at 0x7ff7de0d7e90>
In [5]:
# start simulation
while t<total_time:
t += dt
u = u + dt * (f(u) + D * lap(u,dx) )
# we don't need to plot again, just to update the data of the plot
if len(N) == 1:
p.set_data(x,u)
if len(N) == 2:
p.set_data(u)
ax.set_title("time={:5.1f}".format(t))
fig.canvas.draw()