Code
Chapt 4: PDE
Feng Wan 
B832@仲英楼 or 3133@泓理楼 // School of Physcis @ Xian Jiaotong University
elliptic
Poission eqn
Time-independent Schrodinger
\[ f_{xx} + f_{yy} = 0 \]
\[ \frac{dy}{dx} = \pm i \]
Relaxation methodparabolic
Diffusion eqn
Time-dependent Schrodinger
\[ f_t - \alpha f_{xx} = 0 \]
\[ \frac{dy(x)}{dx(t)} = \infty \]
Code
plt.figure(figsize = (3, 3))
x = np.array([0, 1])
for i in range(6):
y = np.array([1, 1]) * i * 0.1
plt.plot(x, y, 'k-', alpha = 0.3)
for i in range(11):
y = np.array([1, 1]) * i * 0.1
x = np.array([0.5, 1])
plt.plot(y, x, 'k-', alpha = 0.3)
plt.plot([0.5], [0.5], 'ro', lw = 8)
plt.text(0.55, 0.55, r'$P$', fontsize = 15)
plt.axis('off')
plt.show()
Forwardhyperbolic
wave equation
EM wave
\[ f_{tt} - c^2f_{xx} = 0 \]
\[ \frac{dy(x)}{dx(t)} = \pm c \]
Code
plt.figure(figsize = (3, 3))
x = np.array([0, 1])
for i in range(6):
y = np.array([1, 1]) * i * 0.1
plt.plot(x, y, 'k-', alpha = 0.3)
for i in range(11):
y = np.array([1, 1]) * i * 0.1
x = np.array([0.5, 1])
plt.plot(y, x, 'k-', alpha = 0.3)
plt.plot([0.5], [0.5], 'ro', lw = 8)
plt.text(0.5, 0.6, r'$P$', fontsize = 15)
y1 = np.array([0, 0.5])
y2 = np.array([1, 0.5])
x = np.array([0, 0.5])
plt.fill_between(x, y1, y2, color = 'k')
x = np.array([0.5, 1])
y1 = np.array([0.5, 0])
y2 = np.array([0.5, 1])
plt.fill_between(x, y1, y2, color = 'k')
plt.axis('off')
plt.show()
Forward因变量用它们在自变量的许多离散的格点上之值来描述， 通过适当的离散化，偏微分方程就化为一大组差分方程。
4.1 差分法
用差分近似式替代导数，或用适当近似式替代含有导数的表达式，得到这些近似值满足的代数方程
差分法的基本思想 ：用差分代替微分方程中的微分，将连续的变量离散化。通过网络分割法，将函数定义域分成大量的相邻而不重合的子区域，网络线交点成为节点。节点 —— 正则内点，非正则内点
\(\longrightarrow\)
\[h_i = \tau_i = h\]
mesh slicing
:::{#exm-1} 常微分方程
\[ \left\lbrace \begin{aligned}[l] &y'' - x^2y= e^{-x}, ~~ x\in[0, 1] \\ &y(0) = 0, ~~ y(1) = 1 \end{aligned} \right. \] 取 \(h = (b - a)/n = 1/10, ~~ x_i = ih, ~~ i = 0,1,2,...,10\)
:::
利用二阶中心差分公式 \(y_i'' = \frac{y_{i+1} - 2y_i + y_{i-1}}{h^2}\)
\[ \xrightarrow{\text{difference equation}} \left\lbrace \begin{aligned}[l] &y_{i+1} - (2+h^2x_i^2)y_i + y_{i-1} = h^2e^{-x_i} \\ &y_0 = 0, ~~ y_{10} = 1, ~~ i = 1, 2, ..., 9 \end{aligned} \right. \]
椭圆型方程最简单的形式就是Laplace方程 \(\Delta u = u_{xx} + u_{yy} = 0\) 它是如下的Possion方程的特殊情况 \[ \Delta u = u_{xx} + u_{yy} = S(x, y) \] 例如：我们将考虑二维空间 \((x, y)\) 内的场量\(\phi\)的椭圆型边值问题，方程为 \[ -[\partial_x^2 + \partial_y^2]\phi = S(x, y) \]
边界条件：
- Dirichlet问题 \(\phi|_G = g(p)\)
第一类边值问题 - Neumann问题 \(\frac{\partial \phi}{\partial n}\bigg|_G = g(p)\)
第二类边值问题 - 混合问题 \(\phi|_G + g_1(p) \frac{\partial \phi}{\partial n}\bigg|_G = g_2(p)\)
第三类边值问题
取 Dirichlet 型边界条件（第一类边值问题），为方便起见边界线取为单位正方形，边值问题就是要用方程求出单位正方形内每处的\(\phi\)。
首先定义一个网格，覆盖\((x, y)\)平面内的单位正方形。为方便起见，我们取格子间隔 \(h\) (步长)是均匀的，并且在两个方向上相等，使得单位正方形被 \((N+1) \times (N+1)\) 个格点覆盖。这些格点用指标 \((i, j)\) 编号，其中 \(i, j = 0,1,2,3, ..., N\)。
Code
plt.figure(figsize = (5, 5))
x = np.array([0, 1])
for i in range(11):
y = np.array([1, 1]) * i * 0.1
plt.plot(x, y, 'k-', alpha = 0.3)
plt.plot(y, x, 'k-', alpha = 0.3)
plt.text(-0.1, -0.05, r'$x = 0$', fontsize = 15)
plt.text(0.9, -0.05, r'$x = 1$', fontsize = 15)
plt.text(-0.3, -0.05, r'$y = 0$', fontsize = 15)
plt.text(-0.15, 0.97, r'$y = 1$', fontsize = 15)
plt.axis('off')
plt.show()
5.1 将场域离散化
把偏导数用节点上的函数值表示出来，分别以\(f_0, f_1, f_2, f_3, f_4\), 代表在\(0, 1-4\) 的函数值
由1，3点值给出 \(\left( \frac{\partial \phi}{\partial x} \right)_0 \approx \frac{\phi_1 - \phi_0}{h_1} ~~~~~ \left( \frac{\partial \phi}{\partial x} \right)_0 \approx \frac{\phi_0 - \phi_3}{h_3}\)
显然，这种单侧差商误差较大，引入\(a, b\)，由\(f_1\)和\(f_3\) 的泰勒展开式，可构造出下面的关系式（考虑到\(h_1\)与\(h_3\)可能不同）
\[ \alpha(\phi_1 - \phi_0) + \beta(\phi_3 - \phi_0) = \left( \frac{\partial \phi}{\partial x} \right)_0 (\alpha h_1 - \beta h_3) + \frac{1}{2}\left( \frac{\partial^2 \phi}{\partial x^2} \right)_0(\alpha h_1^2 + \beta h_3^2) + ... \]
令第二项系数为零，得
\[ - \left[ \frac{\phi_1 -2\phi_0 + \phi_3}{h^2} + \frac{\phi_2 -2\phi_0 + \phi_4}{h^2} \right] = S_0 \]
\[ \Longrightarrow - \left[ \frac{\phi_{i+1, j} -2\phi_{i,j} + \phi_{i-1,j}}{h^2} + \frac{\phi_{i,j+1} -2\phi_{i,j} + \phi_{i,j-1}}{h^2} \right] = S_{i,j} \]
\[ - \left[ \frac{\phi_{i+1, j} -2\phi_{i,j} + \phi_{i-1,j}}{h^2} + \frac{\phi_{i,j+1} -2\phi_{i,j} + \phi_{i,j-1}}{h^2} \right] = S_{i,j} \]
\[ -[(\delta_i^2\phi)_{ij} + (\delta^2_j\phi)_{ij}] = h^2S_{ij} ~~~~ \Longleftrightarrow ~~~~\left\lbrace \begin{aligned}[l] (\delta_i^2\phi)_{ij} \equiv \phi_{i+1, j} + \phi_{i-1, j} -2\phi_{ij} \\ (\delta_j^2\phi)_{ij} \equiv \phi_{i, j+1} + \phi_{i, j-1} -2\phi_{ij} \end{aligned} \right. \]
5.2 二维问题的差分格式
\[ \left\lbrace \begin{aligned} &u_{xx} + u_{yy} = 0, ~~ (x, y) \in D \\ &u(x, y) = \log[(1 + x)^2 + y^2], ~~ (x, y) \in \partial D \end{aligned} \right. \] \(D = \{ (x, y)|0 \leq x, y \leq 1 \}\)
\[ -[\phi_{i+1, j} + \phi_{i-1, j} + \phi_{i, j+1} + \phi_{i, j-1} -4\phi_{i,j}] = 0 \]
\[ -\left[ \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} \right] \phi = 0 \]
取特殊的网格 \(h=1/3\), 如右图，求差分方程的矩阵形式
Code
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
u = np.zeros((101, 101))
u[-1, :] = 10.0 #边界条件,其他三个边界条件包含在u的初始定义里
uold = u.copy()
u[1:-1, 1:-1] = u[1:-1, 1:-1] + 12.5
nn = 0
while np.max(np.max(np.abs(u-uold)))>=0.00000001:
uold=u.copy()
u[1:-1,1:-1]=0.25*(u[2:,1:-1]+u[0:-2,1:-1]+u[1:-1,2:]+u[1:-1,0:-2])
nn=nn+1
fig = plt.figure(figsize = (10, 4))
ax = fig.add_subplot(121, projection='3d')
x = y = np.linspace(0, 1, 101)
X, Y = np.meshgrid(x, y)
ax.plot_surface(X, Y, u, cmap = 'jet')
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
plt.subplot(122)
plt.imshow(u, cmap = 'jet')
plt.colorbar()
plt.show()
Code
u = np.zeros((301, 201))
x = np.linspace(0, 3, 301)
y = np.linspace(0, 2, 201)
u[-1, :] = np.sin(3.0 * np.pi * y / 2.0)
u[:, -1] = np.sin(3.0 * np.pi * x / 3.0) * np.cos(np.pi * x / 3.0)
uold = u.copy()
u[1:-1, 1:-1] = u[1:-1, 1:-1] + 12.5
nn = 0
while np.max(np.max(np.abs(u-uold)))>=0.00000001:
uold=u.copy()
u[1:-1,1:-1]=0.25*(u[2:,1:-1]+u[0:-2,1:-1]+u[1:-1,2:]+u[1:-1,0:-2])
nn=nn+1
fig = plt.figure(figsize = (10, 4))
ax = fig.add_subplot(121, projection='3d')
X, Y = np.meshgrid(y, x)
ax.plot_surface(X, Y, u, cmap = 'jet')
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
plt.subplot(122)
plt.imshow(u, cmap = 'jet')
plt.colorbar()
plt.show()
Gauss-Seidel 迭代法又可以称为
对于
- \(\omega > 1\) 对应于“超松弛”
- \(\omega < 1\) 则意味“低松弛”
Code
x = np.linspace(0, 1, 100)
y = lambda x: np.exp(-20*(x-0.5)**2)
yy = lambda x: (np.exp(-20*(x-0.5)**2)-np.exp(-20*(x-1.5)**2)-np.exp(-20*(x+0.5)**2))
plt.subplot(121)
plt.plot(x, y(x), 'r-', label = 'y')
plt.plot(x, yy(x), 'b-', label = 'yy')
plt.legend()
plt.subplot(122)
plt.plot(x, y(x) - yy(x))
plt.show()
Code
y_r = lambda x, t: (1.0/np.sqrt(1+80*t)*(np.exp(-20*(x-0.5)**2/(1.0+80*t)) - np.exp(-20*(x-1.5)**2/(1.0+80.0*t))-np.exp(-20*(x+0.5)**2/(1.0+80.0*t))))
t_end = 0.1
t_step = 10000
times = np.linspace(0, t_end, t_step + 1)
dt = times[1] - times[0]
dx = x[1] - x[0]
rr = dt / dx**2.0
yrr = y_r(x, t_end)
y_solve = np.zeros([len(x), t_step + 1])
y_solve[:, 0] = yy(x)
source = np.zeros(x.shape)
def explicit_parabolic(x, y_solve, rr, dt, source):
for nn in range(t_step):
for i in range(1, len(x)-1):
y_solve[i,nn+1]=rr*y_solve[i+1,nn]+(1-2*rr)*y_solve[i,nn]+rr*y_solve[i-1,nn] + dt * source[i]
explicit_parabolic(x, y_solve, rr, dt, source)
for i in range(0, t_step+1, 2000):
plt.plot(x, y_solve[:, i], lw = 5, label = str(i))
plt.plot(x, yrr, 'b-.', lw = 3, label = 'accurate')
plt.legend()
plt.title('r = ' + str(rr))
plt.show()
Code
t_end = 0.1
t_step = 1000
times = np.linspace(0, t_end, t_step + 1)
dt = times[1] - times[0]
dx = x[1] - x[0]
rr = dt / dx**2.0
yrr = y_r(x, t_end)
y_solve = np.zeros([len(x), t_step + 1])
y_solve[:, 0] = yy(x)
source = np.zeros(x.shape)
explicit_parabolic(x, y_solve, rr, dt, source)
for i in range(0, t_step+1, 200):
plt.plot(x, y_solve[:, i], lw = 5, label = str(i))
plt.plot(x, yrr, 'b-.', lw = 3, label = 'accurate')
plt.legend()
plt.title('r = ' + str(rr))
plt.show()
Code
x = np.linspace(0, 20, 500)
t_end = 25.0
t_step = 100000
times = np.linspace(0, t_end, t_step + 1)
dt = times[1] - times[0]
dx = x[1] - x[0]
rr = dt / dx**2.0
print('r = ', rr)
y_solve = np.zeros([len(x), t_step + 1])
for i in range(len(x)):
if x[i] >= 10 and x[i] <= 11:
y_solve[i, 0] = 1.0
source = np.zeros(x.shape)
explicit_parabolic(x, y_solve, rr, dt, source)
fig = plt.figure(figsize = (10, 4))
ax = fig.add_subplot(121, projection='3d')
X, Y = np.meshgrid(times, x)
ax.plot_surface(X, Y, y_solve, cmap = 'jet')
ax.set_xlabel(r'$t$')
ax.set_ylabel(r'$x$')
ax.set_zlabel(r'$\Phi$')
plt.subplot(122)
plt.imshow(y_solve, cmap = 'jet', aspect='auto')
plt.colorbar()
plt.show()r = 0.15562562500000002
Code
x = np.linspace(0, 1, 101)
yy = lambda x: (np.exp(-20*(x-0.5)**2)-np.exp(-20*(x-1.5)**2)-np.exp(-20*(x+0.5)**2))
y_r = lambda x, t: (1.0/np.sqrt(1+80*t)*(np.exp(-20*(x-0.5)**2/(1.0+80*t)) - np.exp(-20*(x-1.5)**2/(1.0+80.0*t))-np.exp(-20*(x+0.5)**2/(1.0+80.0*t))))
t_end = 0.1
t_step = 100
times = np.linspace(0, t_end, t_step + 1)
dt = times[1] - times[0]
dx = x[1] - x[0]
rr = dt / dx**2.0
yrr = y_r(x, t_end)
y_solve = np.zeros([len(x), t_step + 1])
y_solve[:, 0] = yy(x)
def implicit_parabolic(x, y_solve, rr, dt, source):
imat = np.eye(len(x) - 2) * (1 + 2.0 * rr)
for i in range(len(x) - 3):
imat[i, i + 1] = -rr
imat[i + 1, i] = -rr
iimat = np.linalg.inv(imat)
for nn in range(t_step):
temp = y_solve[1:-1, nn] + source[1:-1] * dt
temp[0] = temp[0] + rr*y_solve[0, nn+1]
temp[-1] = temp[-1] + rr*y_solve[-1, nn+1]
y_solve[1:-1, nn + 1] = np.dot(iimat, temp)
source = np.zeros(x.shape)
implicit_parabolic(x, y_solve, rr, dt, source)
for i in range(0, t_step + 1, 20):
plt.plot(x, y_solve[:, i], lw = 5, label = 't = ' + str(times[i]))
plt.plot(x, y_solve[:, -1], lw = 5, color = 'red', label = 'implicit')
plt.plot(x, yrr, 'b-.', lw = 3, label = 'accurate')
plt.legend()
plt.title('r = ' + str(rr))
plt.show()
Code
def avg_implicit_parabolic(x, y_solve, rr, dt, source):
imat = np.eye(len(x) - 2) * (1 + rr)
for i in range(len(x) - 3):
imat[i, i + 1] = -0.5 * rr
imat[i + 1, i] = -0.5 * rr
iimat = np.linalg.inv(imat)
for nn in range(t_step):
temp = 0.5 * rr * y_solve[2:, nn] + (1 - rr) * y_solve[1:-1, nn] + 0.5 * rr * y_solve[0:-2, nn] + source[1:-1] * dt
temp[0] = temp[0] + 0.5 * rr * y_solve[0, nn+1]
temp[-1] = temp[-1] + 0.5 * rr*y_solve[-1, nn+1]
y_solve[1:-1, nn + 1] = np.dot(iimat, temp)
x = np.linspace(0, 1, 101)
yy = lambda x: np.sin(np.pi * x)
y_r = lambda x, t: np.exp(-np.pi**2*t)*np.sin(np.pi*x)
t_end = 0.4
t_step = int(t_end / 0.0025)
times = np.linspace(0, t_end, t_step + 1)
dt = times[1] - times[0]
dx = x[1] - x[0]
rr = dt / dx**2.0
yrr = y_r(x, t_end)
y_solve = np.zeros([len(x), t_step + 1])
y_solve[:, 0] = yy(x)
y_solve[0, 0] = 0.0
y_solve[-1, 0] = 0.0
y_solve2 = y_solve.copy()
y_solve0 = y_solve.copy()
source = np.zeros(x.shape)
avg_implicit_parabolic(x, y_solve2, rr, dt, source)
implicit_parabolic(x, y_solve, rr, dt, source)
explicit_parabolic(x, y_solve0, rr, dt, source)
#plt.plot(x, y_solve0[:, -1], 'r-', label = 'explicit')
plt.plot(x, y_solve[:, -1], 'b-', label = 'implicit')
plt.plot(x, y_solve2[:, -1], 'c-.', label = 'avg')
#plt.plot(x, yrr, 'k', lw = 3, label = 'accurate')
plt.legend()
plt.title('r = ' + str(rr))
plt.show()
显式公式表明，需要两行已知的数据才能求出下一行的数值。 