aix

aix

偏微分方程

  • 学的时候记得笔记都是英语的加上又懒得翻译就只好拖到了现在

General

  • 偏微分方程(PDE, partial differential equation)的解是一个关于时间和空间的方程 $u(x, t)$,我们用 u 代表在特定时间和空间的函数值。
  • 在常微分方程的数值解中,我们是通过将 0-T 划分为若干个 $\Delta t$,通过计算每个 $\Delta t$ 位置的值来计算 $u(t)$。在偏微分方程中,我们构造下面的图形(mesh)来辅助我们的计算:

mesh

我们用 mesh 网格图上的 u 值来近似计算偏微分的值,我们观察临近的矩形,他们构成偏导的近似值。这样我们将偏微分方程的空间离散化,转换为常微分方程的系统。

空间离散化的方案(spatial discretization schemes)有:

  • 有限元法(Finite Element Methods)
  • 有限差分法(Finite Difference Methods)
  • 有限体积法(Finite Volume Methods)

Finite Difference Formulas

我们以一个简单的偏微分方程为例:

k 代表导热系数。

假设网格在空间上的距离为 h,那么我们有:

  • 一阶偏导的前向差分(Forward difference for first spatial derivative)
  • 一阶偏导的向后差分(Backward difference for first spatial derivative)
  • 一阶偏导的中心差分(Central difference for first spatial derivative)
  • 二阶偏导的中心差分(Central difference for second spatial derivative)

Discretization Stencil

  • 我们在网格来近似计算空间函数$u(x, t)$偏导(前向差分,中心差分等等)的这种方法叫做stencil,如果用空间内更多更远的点来近似一个点的偏导值,那么一般情况下误差就会变的更小。但对于有些情况却不同,补充的来说,用更多的点近似来减少误差只对光滑的函数有效,对于那些不连续的函数,相邻点之间的关系很小,使用太多的点反而会影响准确度。
  • 我们用 stencils size 来描述有多少个点被用于近似偏导值。我们的误差和使用的近似方法的阶数成正比关系。我们假设误差为 e,h 为 网格缩小比例,q 是使用方法的阶数,我们有这样的关系:$e \sim h^q$
  • 需要注意:减小空间网格的大小的同时也要减小每一步的时间$\Delta t$来确保稳定性。

Solving PDE

  • 解决 PDE 问题就是解决一个初始边界值问题(initial boundary value problem):$u(x, t)$ 在给定的初始条件和边界条件下是如何变化的。

我们仍以热传导方程$\frac{\partial u}{\partial t} = k \frac{\partial^2 u}{\partial x^2}$为例。在网格内部的点上,我们可以将求偏导转换为计算一组离散化的点,运用我们上面的二阶中心差分的近似计算:

我们想象 i 从 0 到 N-1,总共有 N-1 个类似上面的式子,即我们把 t 时间的 x 范围划分为 N-1 块,每一块都可以通过它临近的点来计算。下面把这 N-1 个式子合成在一起:

这样上面的方程就可以写成类似下面的形式:

f包含边界条件的信息。A 是一个 $N-1 \times N-1$ 的系数矩阵,当x=1和x=N-1的时候都无法计算,所以设置行列式为全 0,需要 f 来辅助确定值以保证边界条件:

这样我们就可以用 ODE 的方法来解决 PDE 问题,比如说我们运用前向欧拉法:

We discretize the spatial domain into N spaces with N+1 grid points. How many unknown variables are evolved in the resulting system of ODEs?

答案是 N-1,$u_0$ 和 $u_{N+1}$ 都是边界点(已知),这样在网格的内部就只有 N-1 个未知的点。

下面我们考虑计算的误差:PDE 方程的解的精确度分为空间和时间上的。我们取 e 代表误差,$\Delta t $ 代表每一步选择迭代的时间差(time_step),h 代表相邻 x 之间的距离,也可以看做网格大小(mesh_size),那么有关系:$e \approx (\Delta t)^p+h^q$。其中 p 和 q 分别是我们选择的 ODE 方法的阶数和偏导数值解方法的阶数。如我们选择 FE (forward eluer) + 中心差分(Central Difference),那么 p 值为 1,q 值为 2。因为它们分别是 1 阶和 2 阶精确度的方法。如果 h 减小一半,那么误差就会变为四分之一。

需要注意,我们前面也提到了,一个小的 h 需要同时一个更小的 $\Delta t$,否则求解将会变得不再稳定。怎么理解呢,我们仍用热传导方程来理解:当我们观察更小的一段的变化的时候,不同段直接的交互会更加的频繁,我们需要用更小的时间来确保可以更加精确的捕获这些变化。

Boundary Conditions

  • Dirichlet conditions

    Dirichlet conditions fix the value of the solution u(x,t) at the boundary. This is like pinning the solution to a certain value at the boundary.

即在边界指定函数的分布形式

  • Neumann conditions

    Neumann conditions fix the value, not of the solution, but of its spatial derivative $\frac{\partial u}{\partial x}$ at the boundary. This is like fixing the stresses at the boundary.

在边界指定外法线方向上的导数的数值

  • Robin conditions

    Robin conditions fix the value of a linear combination of the solution itself and its partial derivative at the boundary.

可以看做是第1和2条件的组合,要求偏导和函数本身的数值。

Linear System

对于转换为 ODE 的 PDE 问题,我们也可以用一些隐式的方法求解,如隐式欧拉法

求解这个方法,移项:

它的每一步相当于求解一个线性系统(linear system)(这不就是线性方程组吗),$M \boldsymbol u=b$
这么做的好处在于,这样让我们可以使用更大的 $\Delta t$ 相比较于显式方法(隐式方法的稳定性,参考 ODE 那节)

求解线性系统有很多的现成方法:

  • A 是一般矩阵(general):高斯消元(Gaussian elimination),$O(N^3)$
  • A 是三角矩阵(triangular):Forward substitution/backward substitution,$O(N^2)$
  • A 是三对角线矩阵(tri-diagonal):托马斯算法(Thomas algorithm),$O(N)$

Iterative Algorithms for Linear Systems

当线性系统中的 M 是一个很大的矩阵,又不特殊时,用高斯消元太慢了并且太占用空间。这时就需要一个更高效的方法:迭代法再次上线。我们先猜一个初始值 v,通过计算 Mv 和 b 比较,修改 v 来一步步逼近真实解。

可选的迭代方法有很多,如 Krylov Subspace Method(krylov 子空间算法),conjugate gradients(共轭梯度法),GMRes(广义最小残差法)。

学一下共轭梯度法然后补充在这里

Nonlinear Systems

对于如下非线性形式的 ODEs:

使用隐式欧拉法得到迭代公式:

求解$u_{n+1}$等同于在这样一个方程中寻找根:$R(u_{n+1})=u_{n+1}-\Delta tf(u_{n+1})=0$,常用的计算方法有:二分法,牛顿法

Sample

我们考虑污染物在一维空间的传播问题。令 $u(x, t)$ 表示在 t 时间 x 位置处污染物的浓度。u 的变化用对流扩散方程来定义:

扩散速率$k = 5 \times 10^{-4} m^2/s$,对流速度$ c= 0.5m/s$
当 $t=0$ 时污染物的分布图像为:

init_ux.png

边界条件是第一类边界条件,在边界的时候$u(x,t)$始终等于 0。

我们采用有限差分法来计算污染物在 $[0, 1]$ 区间内随着时间从 0 变到 1 的变化。首先将偏导全部拆分近似为差分的形式:

然后用上面的方法转化为求 ODE 问题,代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
import numpy as np
import math
import matplotlib.pyplot as plt


def u0(x):
# 高斯函数
mean = 0.2
sigma = 0.025
return (1/160)*np.exp(-1 * ((x - mean) ** 2) / (2 * (sigma ** 2))) / (math.sqrt(2 * np.pi) * sigma)


# init_ux():
# xs = np.linspace(0, 1, 1000)
# ys = [u0(x) for x in xs]
# plt.plot(xs, ys, "b-")
# plt.ylabel(r"$u_0(x)$")
# plt.yticks([round(x, 2) for x in np.linspace(0, 0.1, 10)])
# plt.savefig("init_ux.png")
# plt.show()
sample_num = 5
nSpace = 600
nTime = 1000
xs = np.linspace(0, 1, nSpace-1)
k = 5e-4
c = -0.5
dt = 1/nTime
h = 1/nSpace
A = np.zeros((nSpace-1, nSpace-1))
u = np.zeros((nSpace-1, 1))

for i in range(nSpace-1):
u[i][0] = u0(i*h)

for i in range(1, nSpace-1-1):
A[i][i-1] = k / (h * h) - c / h
A[i][i] = (-2 * k) / (h * h) + c / h
A[i][i+1] = k / (h * h)
# 采用欧拉方法
times = 0
t = 0
u_record = [u.T.tolist()[0]]
while t <= 1:
u = u + dt * A.dot(u)
t += dt
times += 1
# 每 200 次迭代记录一次
if times == nTime/sample_num:
u_record.append(u.T.tolist()[0])
times = 0


# 采用二阶龙格库塔法
u = np.zeros((nSpace-1, 1))
for i in range(nSpace-1):
u[i][0] = u0(i*h)
t = 0
times = 0

while t <= 1:
k1 = A.dot(u)
k2 = A.dot(u+k1*dt)
u = u + 0.5 * dt * (k1 + k2)
t += dt
times += 1
if times == nTime/sample_num:
u_record.append(u.T.tolist()[0])
times = 0

# 绘图
print(len(u_record))
values_fe = [int(i*250/10) for i in range(10)]
values_rk2 = [int(i*250/10) for i in range(10)]
colors_fe = ["#%02x%02x%02x" % (120, int(g), 200)for g in values_fe]
colors_rk2 = ["#%02x%02x%02x" % (200, int(g), 40)for g in values_rk2]

for i in range(sample_num):
if i == 0:
plt.plot(xs, u_record[i], color="#000080", linewidth=2, label="Init u")
elif i == 1:
plt.plot(xs, u_record[i], color=colors_fe[i], linewidth=2, label="FE Method")
else:
plt.plot(xs, u_record[i], color=colors_fe[i], linewidth=2)
for i in range(sample_num):
if i == 0:
plt.plot(xs, u_record[sample_num + i], color=colors_rk2[i], linewidth=2, label="RK2 Method")
else:
plt.plot(xs, u_record[sample_num + i], color=colors_rk2[i], linewidth=2)

plt.title("Evolution of u")
plt.xlabel("x")
plt.ylabel("u(x, t)")
plt.legend(loc="best")
plt.savefig("ux.png")
plt.show()

向后差分 + 显式欧拉 和 向后差分 + 二阶龙格库塔 方法得到的预测结果,颜色深度的变化代表时间的变化,绘制的图形为每隔 200ms 后在 x 上污染物的分布情况。:

ux_FE.png

FE 误差为 8% 左右,如果我们使用 RK2 代替 FE 方法,误差减小到 0.8%。(用的黑盒函数测的误差,暂时还没复现解析解)

可以尝试修改 nSpace 和 nTime 观察效果。当切换方法为向前差分时,很难让结果稳定,对 nSpace 和 nTime 的要求很苛刻。

Author

Ctwo

Posted on

2021-03-02

Updated on

2021-03-04

Licensed under

Comments