In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['figure.dpi'] = 100
plt.rcParams['figure.figsize'] = [10,5]
plt.rcParams['axes.grid'] = True
import pandas as pd
import pprint as pp
import scipy as sc
import numpy as np
import scipy.linalg as lg # SciPy Linear Algebra Library
import scipy.sparse as sp
import scipy.sparse.linalg as sl
from IPython.display import HTML
from matplotlib import animation, rc
Example: two point value linear problem problem¶
$$y''+y=\cos2t\,,\qquad y(0)=1\,,\qquad y(L)=y_2$$Discretization gives $$y_n''\approx \frac{y_{n-1}-2y_n+y_{n+1}}{h^2}\,,\qquad t_n =hn\,,\qquad h = \frac{L}{N-1}=\text{const}$$
Tridiagoanal matrix with equations elements for $0<n<N-1$
$$D_{n*}=\left(0, 0, \ldots, \frac{1}{h^2}, 1-\frac{2}{h^2}, \frac{1}{h^2}, 0, \ldots, 0, 0\right)\,,\quad b_n=\cos t_n$$plus two more (boundary) equations
$$D_{0, *}=(1, 0, 0,\ldots)\,,\qquad b_0=1$$$$D_{N-1, *}=(\ldots, 0, 0, 1)\,,\qquad b_{N-1}=y_2$$In [2]:
help(lg.solve_banded)
In [3]:
"""
* a01 a12 a23 a34 a45 upper diag
a00 a11 a22 a33 a44 a55 diagonal
a10 a21 a32 a43 a54 * lower diag
"""
def banded_example(L, y2, N=1000):
t = np.linspace(0, L, N)
b = np.cos(2*t)
h = (t[-1]-t[0])/(N-1)
u = np.ones(N)/(h*h)
ab = np.array([u, 1-2*u, u])
""" boundary condition implementation """
ab[1, 0] = 1 # a00=1;
ab[0, 1] = 0 # a01=0;
ab[1, -1] = 1 # a[N-1,N-1] = 1
ab[2, -2] = 0 # a[N-1,N-2] = 0
# print(ab)
b[0] = 1
b[-1] = y2
y = lg.solve_banded((1, 1), ab, b);
plt.plot(t, y)
In [4]:
banded_example(5, -1)
In [5]:
help(sp.spdiags)
In [6]:
def sparse_example(L, y2, N=1000):
t = np.linspace(0, L, N)
b = np.cos(2*t)
h = (t[-1]-t[0])/(N-1)
u = np.ones(N)/(h*h)
ab = np.array([u, 1-2*u, u])
""" boundary condition implementation """
ab[1, 0] = 1 # a00=1;
ab[0, 1] = 0 # a01=0;
ab[1, -1] = 1 # a[N-1,N-1] = 1
ab[2, -2] = 0 # a[N-1,N-2] = 0
# print(ab)
b[0] = 1
b[-1] = y2
#y = lg.solve_banded((1, 1), ab, b);
M = sp.spdiags(ab, [1, 0, -1], N, N, format='csc') # sparse 3-diagonal matrix
y = sl.spsolve(M, b);
plt.plot(t, y)
In [7]:
sparse_example(5, -1)
In [8]:
import scipy.sparse as sp
import scipy.sparse.linalg as sl
def Example_BVP(L, y2=2, N=1000, return_equation=False):
t = np.linspace(0, L, N)
b = np.cos(2*t)
h = (t[-1]-t[0])/(N-1)
u = np.ones(N)/(h*h)
ab = np.array([u, 1-2*u, u])
""" boundary condition implementation """
ab[1, 0] = 1 # a00=1;
ab[0, 1] = 0 # a01=0;
ab[1, -1] = 1 # a[N-1,N-1] = 1
ab[2, -2] = 0 # a[N-1,N-2] = 0
# print(ab)
b[0] = 1
b[-1] = y2
#y = lg.solve_banded((1, 1), ab, b);
M = sp.spdiags(ab, [1, 0, -1], N, N, format='csc') # sparse 3-diagonal matrix
if return_equation==True:
return M, b, h
y = sl.spsolve(M, b);
plt.plot(t, y, label=f'$L={L:.2f}, y2={y2}$')
return t, y
In [9]:
plt.subplot(221); s = Example_BVP(L=2*np.pi, y2=2)
plt.subplot(222); s = Example_BVP(L=np.pi, y2=2)
plt.subplot(223); s = Example_BVP(L=2*np.pi, y2=1)
plt.subplot(224); s = Example_BVP(L=np.pi, y2=-1)
In [10]:
for y2 in np.linspace(-1, 3, 9):
s = Example_BVP(L=3, N=1000, y2=y2)
plt.legend()
Out[10]:
General solution has a form $$y(t) = A\cos t+B\sin t-\frac{1}{3}\cos 2t$$ Using $y(0)=1$ we find $A=\frac{4}{3}$ and for initial value problem we have $B=y'(0)$
In [11]:
def solution(t, B):
return 4/3*np.cos(t)+B*np.sin(t)-1/3*np.cos(2*t)
t = np.linspace(0, 8, 1000)
for B in np.linspace(-5, 5, 11):
y = solution(t, B)
s = plt.plot(t, y, label=f'B={B:.2f}')
t = np.array([np.pi, 2*np.pi])
y = solution(t, B)
plt.plot(t, y, 'o', color='black')
plt.plot(t, [2, 2], 'o', color='black', fillstyle='none')
plt.legend(loc='upper right')
Out[11]:
In [12]:
import numpy.linalg as npl
N = 1000
M, b, h = Example_BVP(0.8, 1, return_equation=True, N=N)
rank = npl.matrix_rank(M.todense(), tol=1e-6)
print("Rank for L=0.8:\t", rank)
M, b, h = Example_BVP(np.pi, 1, return_equation=True, N=N)
rank = npl.matrix_rank(M.todense(), tol=1e-6)
print("Rank for L=π:\t", rank)
The rank of the matrix is lower than its size.
When solving $$Mx=b$$ makes some problems it may be better sometimes to solve an equation $$M^TMx=M^Tb$$ which minimizes the following quadratic form $$f(x) = \frac{1}{2}(Mx-b)^T(Mx-b)=\frac{1}{2}x^TM^TMx-x^TM^Tb+b^Tb$$
$$\frac{\partial f(x)}{\partial x}=M^TMx-M^Tb$$In some sense this is a least square method.
For more in-depth analysis use Singular Value Decomposition (SVD)
In [13]:
def pseudo_inverse_example(L, y2):
plt.figure(figsize=(10,3))
plt.subplot(121)
t, y1 = Example_BVP(L, y2)
plt.legend()
M, b, h = Example_BVP(L, y2, return_equation=True)
y = sl.spsolve(M.transpose()*M, M.transpose()*b);
plt.subplot(122)
plt.plot(t, y, c='C2', label = f'pseudo: $L={L:.2f}, y_2={y2}$')
plt.legend()
In [14]:
pseudo_inverse_example(np.pi, -1)
pseudo_inverse_example(np.pi, 1)
In [15]:
N = 20
M, b, h = Example_BVP(2*np.pi, -1, return_equation=True, N=N)
denseM = M.todense()
plt.matshow(denseM)
Out[15]:
In [16]:
plt.matshow((M.transpose()*M).todense())
Out[16]:
Partial differential equations¶
Most often used in physics:
- Wave-type hyperbolic equations $$\partial_{tt}u-\partial_{xx}u +F(u, t, x) = 0 $$
- Parabolic equations (Schroedinger, diffusion) $$i\partial_{t}u+\partial_{xx}u +F(u, t, x) = 0 $$ $$\partial_{t}u-\partial_{xx}u +F(u, t, x) = 0 $$
- Elliptic (Laplace, Poisson, Helmholtz,...) $$\partial_{tt}u+\partial_{xx}u +F(u, t, x) = 0 $$
Each of these types require individual treatment, and describe different physical situations
Partial differential equations - wave-type equations¶
Wave equations describe propagating (relativistic) waves. Nonlinear versions can admit topological solitons:
- integrable sine-Gordon model $$\partial_{tt}u=\partial_{xx}u-\sin u$$
- nonintegrable $\phi^4$ model: $$\partial_{tt}\phi=\partial_{xx}\phi-2\phi(\phi^2-1)$$
Numerical scheme:>br>
- Discretize the spatial dimention $u_n=u(x_n, t)$ (for example $x_n=nh$)
- write a system of $N$ second order (or $2N$ first order) time differential equations $$\partial_{tt}u_n=D_{nm}u_m-F_n$$
- use some time-stepping method to solve the system for given initial conditions $u(x, 0), \partial_t u(x, 0)$ and boundary conditions, for example: $u(0,t)=f(x)$ and $u(L,t)=g(x)$
In [17]:
from numba import jit, njit
import time
@jit
def φ4_equation(u, dx):
D = np.zeros_like(u)
φ = u[2:-2]
fac = (12*dx*dx);
D[2:-2] = (-30/fac)*u[2:-2] + (16/fac)*(u[1:-3]+u[3:-1]) - (1/fac)*(u[0:-4]+u[4:]);
D[2:-2] -= 2*φ*(φ**2-1)
return D
Try not to use Euler mnethod, better Runge-Kutta or simplectic method (for hamiltonian systems). For example $$\partial_t u_n = v_n$$ $$\partial_t v_n = G(u_n)$$ we can use a second order simplectic method
$$u_{n+1/2} = u_n+\frac12 v_n dt$$$$v_{n+1} = v_n+G\left(u_{n+1/2}\right) dt$$$$u_{n+1} = u_{n+1/2}+\frac12 v_n dt$$In [18]:
@jit
def make_step2(u, v, dx, dt):
u += 0.5*v*dt
v += φ4_equation(u, dx)*dt
u += 0.5*v*dt
# return u, v
Important issue: stability.
For explicit methods time step must be adjusted to the highest eigenvalue of $D_{nm}$ (a saw-type function). $$\lambda_{max}\sim (dx)^{-2}$$
In [19]:
dx=0.1
x = np.linspace(0,10*dx, 11)
plt.plot(x, np.cos(np.pi*x/dx))
Out[19]:
- For wave equations
Nothing can propagete faster than light!
- For parabolic equations however:
therefore it is very difficult to solve the equations with high spatial accuracy using explicit methods.
- For elliptic equations: $dx\sim i dy$, there is no propagation!
In [20]:
N=2000
x = np.linspace(-20, 20, N)
dx = x[1]-x[0]
dt = 0.4*dx
In [21]:
@njit
def example_evolution(x, u, v, Tmax, skip=8):
N = len(x)
dx = x[1]-x[0]
dt = 0.4*dx
nsteps = int(Tmax // dt)
frames = int(nsteps//skip+1)
Map = np.zeros((frames, 800), dtype=np.float64)
for n in range(nsteps):
make_step2(u, v, dx, dt)
if(n%skip==0):
Map[n//skip,:] = u[N//2-400:N//2+400]
return Map, n*dt, x
$\phi^4$ equation has solution escribing moving kinks (topological solitons starting tending to $\phi \to\pm 1$ as $x\to \infty$):
$$\phi_K(x,t)=\pm\tanh\left(\gamma(x-vt)\right)\,\qquad\gamma=(1-v^2)^{-1/2}$$Two such kinks (or rather kink-antikink) going from opposite directions can collide
In [22]:
def colliding_kinks(x, V, a=10):
γ = 1/np.sqrt(1-V**2)
u = np.tanh(γ*(x-a))-np.tanh(γ*(x+a))+1
v = V*(np.cosh(γ*(x-a))**(-2)+np.cosh(γ*(x+a))**(-2))
return u, v
In [23]:
N = 1000
x = np.linspace(-N/40, N/40, N)
dx = x[1]-x[0]
u, v = colliding_kinks(x, 0.18)
# short evolution to compile numba
example_evolution(x, u, v, 0.1)
# real calculations of the compiled function
Map, t, x = example_evolution(x, u, v, 100)
In [24]:
def show_map(Map, t, x=x):
dx = x[1]-x[0]
cmap = ['jet', 'ocean', 'coolwarm', 'RdBu', 'seismic', 'terrain', 'twilight', ][5]
plt.imshow(np.array(Map).T, extent=[0, t, -400*dx, 400*dx], aspect = 12*t/800, cmap=cmap, origin='lower')
plt.xlabel('time $t$')
plt.ylabel('$x$')
plt.colorbar(label='$\phi(x,t)$')
In [25]:
N = 1000
x = np.linspace(-N/40, N/40, N)
u, v = colliding_kinks(x, 0.2)
example_evolution(x, u, v, 0.1)
N = int(1e4)
show_map( *example_evolution(x, u, v, 100) )
In [26]:
V = [0.18, 0.20, 0.222, 0.247, 0.25, 0.272]
plt.figure(figsize=(12, 5),dpi=200)
plt.grid(False)
for n, velocity in enumerate(V):
plt.subplot(2, 3, n+1)
plt.grid(False)
u, v = colliding_kinks(x, velocity)
show_map(*example_evolution(x, u, v, 100))
plt.title(f'$v = {velocity}$')
plt.show()
In [27]:
@njit
def evolution_center(x, u, v, Tmax, skip=8):
N = len(x)
dx = x[1]-x[0]
dt = 0.4*dx
nsteps = int(Tmax // dt)
frames = int(nsteps//skip+1)
fi = np.zeros(frames, dtype=np.float64)
t = np.linspace(0, frames*skip*dt, frames)
for n in range(nsteps):
make_step2(u, v, dx, dt)
if(n%skip==0):
fi[n//skip] = u[N//2]
return t, fi
In [28]:
velocity = 0.247
u, v = colliding_kinks(x, velocity)
t, fi = evolution_center(x, u, v, 200)
plt.xlabel("time $t$"); plt.ylabel("field at center $\phi(0,t)$")
a = plt.plot(t, fi)
In [29]:
def make_scan(velocities):
scan = []
for n, vel in enumerate(velocities):
u, v = colliding_kinks(x, velocities[n])
t, fi = evolution_center(x, u, v, 150)
scan.append(fi)
return velocities, t, np.array(scan).T
In [30]:
%%time
N = 3000
x = np.linspace(-150, 150, N)
velocities = np.linspace(0.18, 0.28, 500)
vel1, t1, scan1 = make_scan(velocities, )
In [31]:
%%time
N = 3000
x = np.linspace(-150, 150, N)
velocities = np.linspace(0.224, 0.230, 500)
vel2, t2, scan2 = make_scan(velocities)
In [32]:
def show_scans():
cmap = ['jet', 'ocean', 'coolwarm', 'RdBu', 'seismic', 'terrain', 'twilight', ][5]
fig = plt.subplot(2,1,1)
vel, t, scan = vel1, t1, scan1
plt.imshow(np.array(scan), extent=[vel[0], vel[-1], 0, t[-1]], origin='lower',
aspect = (vel[-1]-vel[0])/100/6, cmap=cmap)
plt.plot([vel2[0], vel2[-1], vel2[-1], vel2[0], vel2[0]],
[0, 0, 149, 149, 0], 'red')
plt.ylabel('time $t$')
plt.xlabel('velocity $v$')
plt.grid(False)
plt.colorbar(label='$\phi$')
fig = plt.subplot(2,1,2)
vel, t, scan = vel2, t2, scan2
plt.imshow(np.array(scan), extent=[vel[0], vel[-1], 0, t[-1]], origin='lower',
aspect = (vel[-1]-vel[0])/100/6, cmap=cmap)
plt.ylabel('time $t$')
plt.xlabel('velocity $v$')
plt.grid(False)
plt.colorbar(label='$\phi$')
In [33]:
show_scans()
Soliton - wave collision¶
In $\phi^4$ model there is a unique interaction between a wave and a soliton.
In [34]:
def wave(x, t, omega=2.2):
k = np.sqrt(omega**2-4)
return np.cos(k*x+omega*t) * 0.5* (1+np.tanh(x-20+k/omega*t))
N = 8000
x = np.linspace(-N/40, N/40, N)
omega, A = 2.5, 0.4
k = np.sqrt(omega**2-4)
u = np.tanh(x) + A * wave(x, 0, omega)
h = 1e-6
fac = (12*h*h);
v = A * ( wave(x, -2*h, omega) - 8*wave(x, -h, omega) + 8*wave(x, h, omega) - wave(x, -2*h, omega))
In [35]:
show_map(*example_evolution(x, u, v, 200))
a =plt.title("Negative radiation pressure i.e. tractor beam")
Elliptic equations:¶
Example: 2-D Poisson equation:
$$\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}=\rho$$with boundary condition $$ \bigl.u\bigr|_{\partial D} = g$$
The lowest order $\mathcal{O}(h^2)$ disretization has the form
$$\frac{u_{i-1, j}+u_{i+1, j}+u_{i, j-1}+u_{i+1, j+1}-4u_{i,j}}{h^2}=\rho_{i,j}$$We have $N^2$ unknowns! so the appropriate matrix would have $N^2\times N^2$ elements.
Fortunatelly most ofthem would be zeros.
In [36]:
def make_grid1(N):
for i in np.arange(N):
for j in np.arange(N):
plt.annotate(f'({i}, {j})', (j+0.1, i+0.1))
if i==0 or i ==N-1 or j==0 or j==N-1:
plt.plot(i, j, 'go')
else:
plt.plot(i, j, 'bo', fillstyle='none')
plt.plot([2, 3, 4], [3,3,3], 'o-', color='black')
plt.plot([3,3,3], [2, 3, 4], 'o-', color='black')
plt.show()
In [37]:
make_grid1(6)
Let's introduce a large multiindex: $0\leq I_{i,j} = iN+j<N^2$.
The Poisson equation reads now
Boundary conditions on a square $0<k<N$:
$$u_{[0:N]}=u_{bottom}\qquad u_{[0:N(N-1):N]}=u_{left}$$$$u_{[N(N-1):N^2]}=u_{top}\qquad u_{[(N-1):N^2:N]}=u_{right}$$In [38]:
def make_grid2(N):
for i in np.arange(N):
for j in np.arange(N):
plt.annotate(f'{i*N+j}', (j+0.1, i+0.1))
if i==0 or i ==N-1 or j==0 or j==N-1:
plt.plot(i, j, 'go')
else:
plt.plot(i, j, 'bo', fillstyle='none')
plt.plot([2, 3, 4], [3,3,3], 'o-', color='black')
plt.plot([3,3,3], [2, 3, 4], 'o-', color='black')
plt.show()
In [39]:
make_grid2(6)
In [40]:
@jit
def poisson_equation(u_boundary, ρ, h):
"""
u_boundary: np.array([left, right, top, bottom]) - boundary conditions
ρ: NxN source
"""
N = u_boundary.shape[1]
"""Numba has some problemw with that """
row = []
col = []
data= []
b = np.zeros(N*N)
for i in np.int32(np.arange(N)):
for j in np.int32(np.arange(N)):
I = i*N+j
if i==0 or i==N-1 or j==0 or j==N-1:
row.append(I); col.append(I); data.append(1)
if i==0: b[I] = u_boundary[3, j]
if i==N-1: b[I] = u_boundary[2, j]
if j==0: b[I] = u_boundary[0, i]
if j==N-1: b[I] = u_boundary[1, i]
else:
row.append(I); col.append(I); data.append(-4)
row.append(I); col.append(I+N) ; data.append(1)
row.append(I); col.append(I-N) ; data.append(1)
row.append(I); col.append(I+1); data.append(1)
row.append(I); col.append(I-1); data.append(1)
b[I] = h*h*ρ[I]
M = sp.csc_matrix((data, (row, col)), shape=(N*N, N*N))
u = sl.spsolve(M, b).reshape(N, N);
if N<=10: return u, M, b
else: return u
In [41]:
def show_matrix_structure(N = 8):
x = np.linspace(0, 1, N)
u_boundary = np.array( [ np.sin(np.pi*x) , np.zeros_like(x), np.zeros_like(x), np.zeros_like(x)] )
ρ = np.zeros(N*N)
ρ[N*N//2+N//2] = -100
u, M, b = poisson_equation(u_boundary, ρ, x[1]-x[0])
M = M.todense()
plt.matshow(M, cmap='bwr', vmin=-4, vmax=4)
blocks = np.arange(0, N*N+1, N)-0.5
plt.xticks(blocks, labels=np.int16(blocks+0.5))
plt.yticks(blocks, labels=np.int16(blocks+0.5))
plt.title("structure of the matrix in the Poisson equation")
plt.show()
return M
In [68]:
M = show_matrix_structure()
Note that the structure of the matrix is multidiagonal: the main diagonal and closest up and down and two more $\pm N$.
This should allow to use sp.spdiags function without a necessity of using loops.
In [43]:
def show_LU(M):
N = np.int32(np.sqrt((M.shape)[0]))
blocks = np.arange(0, N*N+1, N)-0.5
plt.xticks(blocks, labels=np.int16(blocks+0.5))
plt.yticks(blocks, labels=np.int16(blocks+0.5))
P, L, U = lg.lu(M)
plt.subplot(121)
plt.title('lower triangular')
plt.imshow(np.int32(L!=0), cmap='Greys')
plt.xticks(blocks, labels=np.int16(blocks+0.5))
plt.yticks(blocks, labels=np.int16(blocks+0.5))
plt.subplot(122)
plt.title('upper triangular')
plt.imshow(np.int32(U!=0), cmap='Greys')
plt.xticks(blocks, labels=np.int16(blocks+0.5))
plt.yticks(blocks, labels=np.int16(blocks+0.5))
print("Nonzero elements of M: ", (M[M!=0]).size)
print("Nonzero elements of LU:",(L[L!=0]).size + (U[U!=0]).size)
plt.show()
In [44]:
show_LU(M)
In [45]:
N = 200
x = np.linspace(0, 1, N)
u_boundary = np.array( [ np.sin(np.pi*x) , np.zeros_like(x), np.sin(2*np.pi*x), np.zeros_like(x)] )
ρ = np.zeros(N*N)
ρ[N*N//2+N//2] = -N**2
u = poisson_equation(u_boundary, ρ, x[1]-x[0])
In [46]:
#plt.imshow(u, extent=[0, 1, 0, 1], cmap='RdBu', vmin=-1, vmax=1)
plt.imshow(u, extent=[0, 1, 0, 1], cmap='Spectral')
plt.colorbar()
levels = np.linspace(-1, 1, 21)
contour = plt.contour(x, 1-x, u, levels, colors='k')
plt.grid(False)
In [47]:
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
X, Y = np.meshgrid(x, 1-x)
p =ax.plot_surface(X, Y, u, cmap='Spectral')
Diffusion equation¶
Due to the stability reasons we have to use implicit method
Crank-Nicholson method for diffusion equation¶
$$u_t=F(u, x, t, u_x, u_{xx})=u_{xx}$$We can use Adams-Moulton trapezoidal rule $u(x_n, t_k)\equiv u^k_n$ with $x_{n+1}-x_n=h$, $t_{n+1}-t_n=\lambda h$
$$u^{k+1}_{n}=u^k_n+\frac{\lambda h}{2}\left(F^k_n+F^{k+1}_n\right)$$For the pure diffusion equation in $\mathcal{O}(h^2)$ we have
$$u^{k+1}_{n}=u^k_n+\frac{\lambda}{2h}\left(u^{k}_{n-1}-2u^{k}_n+u^{k}_{n+1} + u^{k+1}_{n-1}-2u^{k+1}_n+u^{k+1}_{n+1}\right)$$The equation for values at $t_{k+1}$ time can be written in terms of matrices:
$$M_{nm}u^{k+1}_m=b_n$$with
$$M_{nm}=\left(1+\frac{\lambda}{h}\right)\delta_{nm}-\frac{\lambda}{2h}\left(\delta_{n, m+1}+\delta_{n, m-1}\right)$$$$b_n=\left(1-\frac{\lambda}{h}\right)u^k_n+\frac{\lambda}{2h}\left(u^k_{n+1}+u^k_{n-1}\right)$$But beware of boundary conditions for $n=0$ and $n=N-1$!
In [69]:
@jit
def Crank_Nicholson_steps(𝜆_ℎ, u, nsteps=10, ret_skip=1):
v = 0.5*𝜆_ℎ*np.ones_like(u)
tridiag_data = np.array([-v, 1+2*v, -v])
""" boundary condition implementation """
tridiag_data[1, 0] = 1
tridiag_data[0, 1] = 0
tridiag_data[1, -1] = 1
tridiag_data[2, -2] = 0
res = [u]
t = [ 0. ]
for n in range(nsteps):
b = np.empty_like(u)
b[1:-1] = (1-𝜆_ℎ)*u[1:-1] + 0.5*𝜆_ℎ*(u[2:]+u[0:-2])
b[0] = u[0]
b[-1] = u[-1]
u = lg.solve_banded((1, 1), tridiag_data, b);
if(n%ret_skip==0):
res.append(u)
t.append(n*1.0)
return t, np.array(res).T
In [80]:
𝜆_ℎ = 20
N = 1000
x = np.linspace(-1, 1, N)
u = -x#np.exp(-16*x**2)
u[u<0] = 0
In [81]:
t, res = Crank_Nicholson_steps(𝜆_ℎ, u, nsteps=1, ret_skip = 1)
In [82]:
t, res = Crank_Nicholson_steps(𝜆_ℎ, u, nsteps=4000, ret_skip = 400)
res = np.array(res)
plt.plot(x, res)
plt.show()
Split step method - nonlinear Schroedinger (Gross–Pitaevskii) equation¶
$$iu_t = -\frac12u_{xx}+\kappa |u^2|u\,,\qquad \kappa=\pm 1$$with knowns soliton solutions $\kappa=1$: $$u(x, t) = \frac{Ae^{i\Omega t}}{\cosh(Ax)}\,,\qquad \Omega = \frac{A^2}{2}$$
Equation written in the form:
$$iu_t=\left(-\frac12\partial_{xx}+\kappa|u^2|\right)u=(D+V)u$$has a formal solution
$$u(t) = \exp\left[\displaystyle -i\int_{t_0}^tD+V\,dt\right]u(t_0)$$For short times
$$u(t+dt) = \exp\left[-i(D+V)dt\right]\,u(t)$$Note that exponantioation is easily done in eigenbasis:
$$\exp A=O^\dagger \exp\left(OAO^\dagger\right) O$$$V$ is already diagonal, but $D\sim\hat p^2$ is diagonal in a momentum space.
To go to momentum space we need Fourier Transform (especially if $u$ is spatially periodic)
Using Baker–Campbell–Hausdorff formula $$e^Xe^Y=e^{X+Y+{\frac {1}{2}}[X,Y]+{\frac {1}{12}}[X,[X,Y]]-{\frac {1}{12}}[Y,[X,Y]]+\cdots \,,}$$
Therefore in $\mathcal{O}(dt)$ we can split
$$u(t+dt) = e^{-iDdt}e^{-iVdt}\,u(t)$$or in $\mathcal{O}(dt^2)$:
$$u(t+dt) = e^{-\frac{i}{2}Vdt}e^{-iDdt}e^{-\frac{i}{2}Vdt}\,u(t)$$Which can be calculated in the following way
$$u(t+dt) = \mathcal{F}^{-1}\left[e^{-i\tilde Ddt}\mathcal{F}\left(e^{-iVdt}\,u(t)\right)\right]$$second order $\mathcal{O}(dt^2)$ scheme can also be applied:
$$u(t+dt) = e^{-\frac{i}{2}Vdt}\mathcal{F}^{-1}\left[e^{-i\tilde Ddt}\mathcal{F}\left(e^{-\frac{i}{2}Vdt}\,u(t)\right)\right]$$Note that $$Du = \mathcal{F}^{-1}\left[-k^2\mathcal{F}(u)\right]$$ $$e^{Du} = \mathcal{F}^{-1}\left[e^{-k^2}\mathcal{F}(u)\right]$$
In [51]:
def test_spectral_derivative():
N = 1000
x = np.linspace(-20, 20, N)
A = 1
u = A/np.cosh(A*x)
k = np.empty(N)
k[0:N//2+1] = np.arange(N//2+1)*2*np.pi/(x[-1]-x[0])*(N-1)/(N)
n = np.arange(1,N//2)
k[-n] = -k[n]
u_fft = -k*k*np.fft.fft(u)
Du = np.fft.ifft(u_fft)
Du_exact = A**2*(np.sinh(A*x)**2-1)/np.cosh(A*x)**3
plt.semilogy(x, np.abs(Du_exact-np.real(Du)))
plt.title("Spectral derivative error")
plt.xlabel("x")
plt.ylabel("$|Du(x)-u'(x)|$")
In [52]:
test_spectral_derivative()
In [98]:
def split_step_evolution(u, x, dt, nsteps, snaps=0):
N = len(u)
k = np.empty(N)
k[0:N//2+1] = np.arange(N//2+1)*2*np.pi/(x[-1]-x[0])*(N-1)/(N)
n = np.arange(1,N//2)
k[-n] = -k[n]
res = []
t = []
if snaps != 0:
profiles =[]
for n in range(nsteps):
#print(u[N//2])
u = np.exp(-1j*np.conj(u)*u*dt)*u
u_fft = np.fft.fft(u)
u = np.fft.ifft(np.exp(0.5j*k*k*dt)*u_fft)
res.append(np.real(u[N//2:N//2+100:20]))
t.append(n*dt)
if snaps!=0 and n%snaps==0:
profiles.append(u)
if snaps==0:
return np.array(t), np.array(res)
else:
return np.array(t), np.array(res), np.array(profiles).T
In [99]:
def nonlinear_schroedinger_example():
"""Test soliton solution $u(x,t) = \frac{A}{\cosh(Ax) e^{iA^2t/2}}$ """
N = 1001
x = np.linspace(-20, 20, N)
A = 1
u = A/np.cosh(A*x)
h = x[1]-x[0]
dt = h
t, res = split_step_evolution(u, x, dt, 1000)
plt.plot(np.array(t)/np.pi, res)
plt.title("Example solution of NLSE")
plt.xlabel("$t/\pi$")
plt.ylabel("Re $u(x,t)$")
labs = [ 'x = {:.2f}'.format(X) for X in x[N//2:N//2+100:20] ]
plt.legend(labs, title='x')
In [100]:
nonlinear_schroedinger_example()
In [101]:
def NLSE_collision():
N = 2000
x = np.linspace(-20, 20, N)
A = [1, 1.5]
u = A[0]/np.cosh(A[0]*(x+10))*np.exp(-0.5j*x) + A[1]/np.cosh(A[1]*(x-5))*np.exp(0.1j*x)
h = x[1]-x[0]
dt = h
t, res, profiles = split_step_evolution(u, x, dt, 2000, 100)
for k in range(profiles.shape[1]):
a = plt.plot(x, np.abs(profiles[:,k])-k)
a = plt.fill_between(x, np.abs(profiles[:,k])-k, -k-1, alpha = 0.25)
plt.xlabel("x")
plt.title("Soliton collisions in NLSE")
plt.ylabel('$u(x,t)+\\alpha t$')
plt.show()
Example 2: collision of two solitons with initial cionitions: $$u(x, t=0) = \frac{A_0e^{-ivx_0}}{\cosh(A_0(x-x_0))} + \frac{A_1e^{-ivx_1}}{\cosh(A_1(x-x_1))}$$
In [102]:
NLSE_collision()
In [103]:
N = 2000
x = np.linspace(-20, 20, N)
A = [1, 1.5]
u = A[0]/np.cosh(A[0]*(x+10))*np.exp(-0.2j*x) + A[1]/np.cosh(A[1]*(x-10))*np.exp(0.3j*x)
h = x[1]-x[0]
dt = h
t, res, profiles = split_step_evolution(u, x, dt, 4000, 10)
In [89]:
from IPython.display import HTML
from matplotlib import animation, rc
fig, ax = plt.subplots()
def animate(i):
ax.clear()
ax.set_xlim(( -20, 20))
ax.set_ylim((-1, 4))
y = np.abs(profiles[:,i])**2
a = ax.fill_between(x, y, 0, color = "green", alpha=0.2)
a = ax.plot(x, y, "black")
# call the animator. blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, animate, frames=400, interval=20, blit=False)
#HTML(anim.to_html5_video())
#anim.save('NLSE.mp4', writer='ffmpeg', fps=10, dpi=100)
In [60]:
HTML('<center><video controls autoplay loop src="NLSE.mp4" width=100%/></center>')
Out[60]:
Tasks¶
- Analyze collisions of kinks in the sine-Gordon model:
$$\partial_{tt}u=\partial_{xx}u-\sin u$$
Static kink(antykink) solytions are given by
$$u(x, t) = 4\textrm{arctan}\, e^{\pm x}$$
Boundary conditions have to fulfill $\cos(u(\pm L)=1$
- Change functions for boundary conditions in the Poisson equation.
- Solve a Poisson equation on a triangle
- Solve anti-diffusion equation $u_t=-u_{xx}$. What is wrond with the equation?
- How does the soliton collision looks like in a nonlinear Schrodinger equation with a higher nonlinearity? $$i\partial_tu = -\frac12\partial_{xx}u+ |u^3|u$$