8  The time-dependent Ginzburg-Landau equation

Finite differences method applied to a beautiful equation

8.1 Introduction

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: - Euler method to advance time - Five-point stencil to compute the laplacian, periodic boundary conditions are assumed. See an example of the output here: https://www.youtube.com/watch?v=JgE9Px7zsQE

8.2 Code

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import time as tm

# 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

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

initialize and start plotting

plt.ion()
fig = plt.figure(1,figsize=(7,6))
plt.clf()
ax = fig.add_subplot(111)
# 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)
    colorbar_ax = divider.append_axes("right", size="5%", pad=0.15)
    cbar = fig.colorbar(p, cax=colorbar_ax, ticks=[-1,-0.5,0,0.5,1])
ax.set_title("time={:5.1f}".format(0.0))

Text(0.5, 1.0, 'time=  0.0')

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()