七、偏微分方程数值解法
偏微分方程一般可分为抛物型方程、双曲型方程和椭圆型方程三种类型,其一般形式如下
F(uxx,uxy,uyy,ux,uy,u,x,y)=Auxx+Buxy+Cuyy+F~(ux,uy,u,x,y)=0△=B2−4AC{=0⇒抛物型>0⇒双曲型<0⇒椭圆型
\begin{aligned}
&F(u_{xx},u_{xy},u_{yy},u_x,u_y,u,x,y)= Au_{xx}+Bu_{xy}+Cu_{yy}+\tilde{F}(u_x,u_y,u,x,y) =0 \\
&\triangle = B^2 - 4AC \left\{
\begin{aligned}
&=0 \Rightarrow 抛物型 \\
&\gt 0 \Rightarrow 双曲型 \\
&\lt0 \Rightarrow 椭圆型
\end{aligned}
\right.
\end{aligned}
F(uxx,uxy,uyy,ux,uy,u,x,y)=Auxx+Buxy+Cuyy+F~(ux,uy,u,x,y)=0△=B2−4AC⎩⎪⎨⎪⎧=0⇒抛物型>0⇒双曲型<0⇒椭圆型
用差分方法求微分方程近似解时常用公式(通过泰勒展开即可得到)
g′(x0)=1h[g(x0+h)−g(x0)]−h2g′′(ξ1), ξ1∈(x0,x0+h)(0.1)g′(x0)=1h[g(x0)−g(x0−h)]+h2g′′(ξ2), ξ2∈(x0−h,x0)(0.2)g′(x0)=1h[g(x0+h2)−g(x0−h2)]−h224g′′′(ξ3), ξ3∈(x0−h2,x0+h2)(0.3)g′′(x0)=1h2[g(x0+h)−2g(x0)+g(x0−h)]−h212g(4)(ξ4), ξ4∈(x0−h,x0+h)(0.4)g(x0)=12[g(x0−h)+g(x0+h)]−h22g′′(ξ5), ξ5∈(x0−h,x0+h)(0.5)
\begin{aligned}
&g^{'}(x_0) = \frac{1}{h}[g(x_0+h)-g(x_0)]-\frac{h}{2}g^{''}(\xi_1),\ \xi_1 \in (x_0, x_0+h) &(0.1)\\
&g^{'}(x_0) = \frac{1}{h}[g(x_0)-g(x_0-h)]+\frac{h}{2}g^{''}(\xi_2),\ \xi_2 \in (x_0-h, x_0) &(0.2)\\
&g^{'}(x_0) = \frac{1}{h}[g(x_0+\frac{h}{2})-g(x_0-\frac{h}{2})]-\frac{h^2}{24}g^{'''}(\xi_3),\ \xi_3 \in(x_0-\frac{h}{2}, x_0+\frac{h}{2}) &(0.3)\\
&g^{''}(x_0) = \frac{1}{h^2}[g(x_0+h)-2g(x_0)+g(x_0-h)]-\frac{h^2}{12}g^{(4)}(\xi_4),\ \xi_4 \in (x_0-h, x_0+h) &(0.4)\\
&g(x_0) = \frac{1}{2}[g(x_0-h)+g(x_0+h)]-\frac{h^2}{2}g^{''}(\xi_5),\ \xi_5 \in (x_0-h, x_0+h) &(0.5)
\end{aligned}
g′(x0)=h1[g(x0+h)−g(x0)]−2hg′′(ξ1), ξ1∈(x0,x0+h)g′(x0)=h1[g(x0)−g(x0−h)]+2hg′′(ξ2), ξ2∈(x0−h,x0)g′(x0)=h1[g(x0+2h)−g(x0−2h)]−24h2g′′′(ξ3), ξ3∈(x0−2h,x0+2h)g′′(x0)=h21[g(x0+h)−2g(x0)+g(x0−h)]−12h2g(4)(ξ4), ξ4∈(x0−h,x0+h)g(x0)=21[g(x0−h)+g(x0+h)]−2h2g′′(ξ5), ξ5∈(x0−h,x0+h)(0.1)(0.2)(0.3)(0.4)(0.5)
抛物型方程的差分方法
考虑下面的定解问题
{∂u∂t−a∂2u∂x2=f(x,t), x∈(0,l), t∈(0,T)(1.1)u(x,0)=φ(x), x∈[0,l](1.2)u(0,t)=α(t), u(l,t)=β(t), t∈[0,T](1.3)
\left\{
\begin{aligned}
&\frac{\partial u}{\partial t} - a \frac{\partial^2 u}{\partial x^2} = f(x,t),\ x \in (0, l),\ t \in (0, T) &(1.1)\\
&u(x,0)=\varphi(x),\ x \in [0, l] &(1.2)\\
&u(0, t) = \alpha(t),\ u(l, t) = \beta(t),\ t \in [0, T] &(1.3)
\end{aligned}
\right.
⎩⎪⎪⎪⎨⎪⎪⎪⎧∂t∂u−a∂x2∂2u=f(x,t), x∈(0,l), t∈(0,T)u(x,0)=φ(x), x∈[0,l]u(0,t)=α(t), u(l,t)=β(t), t∈[0,T](1.1)(1.2)(1.3)
其中 aaa 是正常数,fff、φ\varphiφ、α\alphaα、β\betaβ 是已知函数,满足 φ(0)=α(0)\varphi(0) = \alpha(0)φ(0)=α(0),φ(l)=β(0)\varphi(l) = \beta(0)φ(l)=β(0)。假设上述定解问题存在唯一解 u(x,t)u(x,t)u(x,t),且具有一定的光滑性
网格剖分
记 D={(x,t)∣0≤x≤l,0≤t≤T}D = \{(x,t)|0 \le x \le l,0 \le t \le T \}D={(x,t)∣0≤x≤l,0≤t≤T},h=lMh=\frac{l}{M}h=Ml,τ=TN\tau = \frac{T}{N}τ=NT,xi=ih(0≤i≤M)x_i=ih(0 \le i \le M)xi=ih(0≤i≤M),tk=kτ(0≤k≤N)t_k = k\tau(0 \le k \le N)tk=kτ(0≤k≤N),用两簇平行直线
x=xi, i=0,1,..,Mt=tk, k=0,1,...,N
x = x_i,\ i=0,1,..,M \\
t = t_k,\ k=0,1,...,N
x=xi, i=0,1,..,Mt=tk, k=0,1,...,N
将区域 DDD 分割成矩形网格。hhh 和 τ\tauτ 称为空间步长和时间步长,网格点 (xi,tk)(x_i, t_k)(xi,tk) 称为节点。记 D‾h={(xi,tk)∣0≤i≤M,0≤k≤N}\overline{D}_h = \{(x_i, t_k)|0 \le i \le M, 0 \le k \le N \}Dh={(xi,tk)∣0≤i≤M,0≤k≤N},Dh={(xi,tk)∣1≤i≤M−1,1≤k≤N}D_h = \{(x_i, t_k)| 1 \le i \le M-1, 1 \le k \le N \}Dh={(xi,tk)∣1≤i≤M−1,1≤k≤N},Γh=D‾h\Gamma_h = \overline{D}_hΓh=Dh\ DhD_hDh。称 DhD_hDh 为内部节点,Γh\Gamma_hΓh 为边界节点。在位于 t=tkt=t_kt=tk 上的所有节点称为第 kkk 层节点。u(xi,tk)u(x_i, t_k)u(xi,tk) 近似值可表示为 uiku_i^kuik
考虑在点 (xi,tk)(x_i, t_k)(xi,tk) 处的方程
∂u∂t(xi,tk)−a∂2u∂x2(xi,tk)=f(xi,tk)(1.4)
\frac{\partial u}{\partial t}(x_i, t_k) - a\frac{\partial^2u}{\partial x^2}(x_i, t_k) = f(x_i, t_k)\quad (1.4)
∂t∂u(xi,tk)−a∂x2∂2u(xi,tk)=f(xi,tk)(1.4)
- 利用 (0.1)(0.1)(0.1)、(0.4)(0.4)(0.4)、初始条件 (1.2)(1.2)(1.2) 和边界条件 (1.3)(1.3)(1.3) 可得古典显格式(向前 Euler 格式)
- 利用 (0.2)(0.2)(0.2)、(0.4)(0.4)(0.4)、初始条件 (1.2)(1.2)(1.2) 和边界条件 (1.3)(1.3)(1.3) 可得古典隐格式(向后 Euler 格式)
- 利用 (0.3)(0.3)(0.3)、(0.4)(0.4)(0.4)、初始条件 (1.2)(1.2)(1.2) 和边界条件 (1.3)(1.3)(1.3) 可得 Richardson 格式
古典显格式
{1τ(uik+1−uik)−ah2(ui+1k−2uik+ui−1k)=f(xi,tk), 1≤i≤M−1,0≤k≤N−1(1.6)ui0=φ(xi), 1≤i≤M−1(1.7)u0k=α(tk), uMk=β(tk), 0≤k≤N(1.8) \left\{ \begin{aligned} &\frac{1}{\tau}(u_i^{k+1}-u_i^k)-\frac{a}{h^2}(u_{i+1}^k-2u_i^k+u_{i-1}^k) = f(x_i, t_k), \ 1 \le i \le M-1, 0 \le k \le N-1 &(1.6)\\ &u_i^0 = \varphi(x_i),\ 1 \le i \le M-1 &(1.7) \\ &u_0^k = \alpha(t_k),\ u_M^k = \beta(t_k),\ 0 \le k \le N &(1.8) \end{aligned} \right. ⎩⎪⎪⎪⎨⎪⎪⎪⎧τ1(uik+1−uik)−h2a(ui+1k−2uik+ui−1k)=f(xi,tk), 1≤i≤M−1,0≤k≤N−1ui0=φ(xi), 1≤i≤M−1u0k=α(tk), uMk=β(tk), 0≤k≤N(1.6)(1.7)(1.8)
记 r=aτh2r = \frac{a\tau}{h^2}r=h2aτ,称 rrr 是步长比,则 (1.6)(1.6)(1.6) 式可改写为
uik+1=(1−2r)uik+r(ui−1k+ui+1k)+τf(xi,tk)
u_i^{k+1} = (1-2r)u_i^k + r(u_{i-1}^k+u_{i+1}^k) + \tau f(x_i, t_k)
uik+1=(1−2r)uik+r(ui−1k+ui+1k)+τf(xi,tk)
从而差分格式 (1.6)−(1.8)(1.6)-(1.8)(1.6)−(1.8) 可改写成向量和矩阵形式
[u1k+1u2k+1⋮uM−2k+1uM−1k+1]=[1−2rrr1−2rr⋱⋱⋱r1−2rrr1−2r][u1ku2k⋮uM−2kuM−1k]+[τf(x1,tk)+rα(tk)τf(x2,tk)⋮τf(xM−2,tk)τf(xM−1,tk)+rβ(tk)]k=0,1,...,N−1
\left[
\begin{matrix}
u_1^{k+1} \\
u_2^{k+1} \\
\vdots \\
u_{M-2}^{k+1} \\
u_{M-1}^{k+1}
\end{matrix}
\right] =
\left[
\begin{matrix}
1-2r & r \\
r & 1-2r &r \\
&\ddots & \ddots & \ddots \\
&&r &1-2r & r \\
&&&r & 1-2r
\end{matrix}
\right]
\left[
\begin{matrix}
u_1^k \\
u_2^k \\
\vdots \\
u_{M-2}^k \\
u_{M-1}^k
\end{matrix}
\right] +
\left[
\begin{matrix}
\tau f(x_1, t_k)+r\alpha(t_k) \\
\tau f(x_2, t_k) \\
\vdots \\
\tau f(x_{M-2}, t_k) \\
\tau f(x_{M-1}, t_k) + r\beta(t_k)
\end{matrix}
\right]
\quad k = 0,1,...,N-1
⎣⎢⎢⎢⎢⎢⎡u1k+1u2k+1⋮uM−2k+1uM−1k+1⎦⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡1−2rrr1−2r⋱r⋱r⋱1−2rrr1−2r⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎡u1ku2k⋮uM−2kuM−1k⎦⎥⎥⎥⎥⎥⎤+⎣⎢⎢⎢⎢⎢⎡τf(x1,tk)+rα(tk)τf(x2,tk)⋮τf(xM−2,tk)τf(xM−1,tk)+rβ(tk)⎦⎥⎥⎥⎥⎥⎤k=0,1,...,N−1
其截断误差如下
Rik(1)=τ2∂2u∂t2(xi,ηik)−ah212∂4u∂x4(ξik,tk), ηik∈(tk,tk+1), ξik∈(xi−1,xi+1)
R_{ik}^{(1)} = \frac{\tau}{2}\frac{\partial^2 u}{\partial t^2}(x_i, \eta_i^k) - \frac{ah^2}{12}\frac{\partial^4 u}{\partial x^4}(\xi_i^k, t_k),\ \eta_i^k \in (t_k, t_{k+1}),\ \xi_i^k \in (x_{i-1}, x_{i+1})
Rik(1)=2τ∂t2∂2u(xi,ηik)−12ah2∂x4∂4u(ξik,tk), ηik∈(tk,tk+1), ξik∈(xi−1,xi+1)
古典隐格式
{1τ(uik−uik−1)−ah2(ui+1k−2uik+ui−1k)=f(xi,tk), 1≤i≤M−1,1≤k≤N(1.11)ui0=φ(xi), 1≤i≤M−1(1.12)u0k=α(tk), uMk=β(tk), 0≤k≤N(1.13) \left\{ \begin{aligned} &\frac{1}{\tau}(u_i^{k}-u_i^{k-1})-\frac{a}{h^2}(u_{i+1}^k-2u_i^k+u_{i-1}^k) = f(x_i, t_k), \ 1 \le i \le M-1, 1 \le k \le N &(1.11)\\ &u_i^0 = \varphi(x_i),\ 1 \le i \le M-1 &(1.12) \\ &u_0^k = \alpha(t_k),\ u_M^k = \beta(t_k),\ 0 \le k \le N &(1.13) \end{aligned} \right. ⎩⎪⎪⎪⎨⎪⎪⎪⎧τ1(uik−uik−1)−h2a(ui+1k−2uik+ui−1k)=f(xi,tk), 1≤i≤M−1,1≤k≤Nui0=φ(xi), 1≤i≤M−1u0k=α(tk), uMk=β(tk), 0≤k≤N(1.11)(1.12)(1.13)
将其改写成矩阵形式
[1+2r−r−r1+2r−r⋱⋱⋱−r1+2r−r−r1+2r][u1ku2k⋮uM−2kuM−1k]=[u1k−1u2k−1⋮uM−2k−1uM−1k−1]+[τf(x1,tk)+rα(tk)τf(x2,tk)⋮τf(xM−2,tk)τf(xM−1,tk)+rβ(tk)]k=1,2,...,N
\left[
\begin{matrix}
1+2r & -r \\
-r & 1+2r &-r \\
&\ddots &\ddots &\ddots \\
&&-r &1+2r & -r \\
&&&-r & 1+2r
\end{matrix}
\right]\left[
\begin{matrix}
u_1^{k} \\
u_2^{k} \\
\vdots \\
u_{M-2}^{k} \\
u_{M-1}^{k}
\end{matrix}
\right] =
\left[
\begin{matrix}
u_1^{k-1} \\
u_2^{k-1} \\
\vdots \\
u_{M-2}^{k-1} \\
u_{M-1}^{k-1}
\end{matrix}
\right] +
\left[
\begin{matrix}
\tau f(x_1, t_k)+r\alpha(t_k) \\
\tau f(x_2, t_k) \\
\vdots \\
\tau f(x_{M-2}, t_k) \\
\tau f(x_{M-1}, t_k) + r\beta(t_k)
\end{matrix}
\right]
\quad k = 1,2,...,N
⎣⎢⎢⎢⎢⎡1+2r−r−r1+2r⋱−r⋱−r⋱1+2r−r−r1+2r⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎡u1ku2k⋮uM−2kuM−1k⎦⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎡u1k−1u2k−1⋮uM−2k−1uM−1k−1⎦⎥⎥⎥⎥⎥⎤+⎣⎢⎢⎢⎢⎢⎡τf(x1,tk)+rα(tk)τf(x2,tk)⋮τf(xM−2,tk)τf(xM−1,tk)+rβ(tk)⎦⎥⎥⎥⎥⎥⎤k=1,2,...,N
其截断误差如下
Rik(2)=−τ2∂2u∂t2(xi,ηik)−ah212∂4u∂x4(ξik,tk), ηik∈(tk−1,tk), ξik∈(xi−1,xi+1)
R_{ik}^{(2)} = -\frac{\tau}{2}\frac{\partial^2 u}{\partial t^2}(x_i, \eta_i^k) - \frac{ah^2}{12}\frac{\partial^4 u}{\partial x^4}(\xi_i^k, t_k),\ \eta_i^k \in (t_{k-1}, t_k),\ \xi_i^k \in (x_{i-1}, x_{i+1})
Rik(2)=−2τ∂t2∂2u(xi,ηik)−12ah2∂x4∂4u(ξik,tk), ηik∈(tk−1,tk), ξik∈(xi−1,xi+1)
Richardson
- Richardson 格式是一个 3 层格式,需要知道 (k-1) 层和 k 层的值才能求出第 (k+1) 层的值
- Richardson 格式是完全不稳定的格式,无实用价值
Crank-Nicolson
为了提高时间方向的精度,记 tk+τ2=tk+τ2t_{k+\frac{\tau}{2}}=t_k + \frac{\tau}{2}tk+2τ=tk+2τ,考虑点 (xi,tk+τ2)(x_i, t_{k+\frac{\tau}{2}})(xi,tk+2τ) 处的微分方程
∂u∂t(xi,tk+τ2)−a∂2u∂x2(xi,tk+τ2)=f(xi,tk+τ2)
\frac{\partial u}{\partial t}(x_i, t_{k+\frac{\tau}{2}}) - a\frac{\partial^2u}{\partial x^2}(x_i, t_{k+\frac{\tau}{2}}) = f(x_i, t_{k+\frac{\tau}{2}})
∂t∂u(xi,tk+2τ)−a∂x2∂2u(xi,tk+2τ)=f(xi,tk+2τ)
利用中心差商 (0.3)(0.3)(0.3)、平均公式 (0.5)(0.5)(0.5) 和二阶差商公式 (0.4)(0.4)(0.4) 可得 Crank-Nicolson 格式
{uik+1−uikτ−a⋅12[ui+1k−2uik+ui−1kh2+ui+1k+1−2uik+1+ui−1k+1h2]=f(xi,tk+τ2), 1≤i≤M−1,0≤k≤N−1(1.19)ui0=φ(xi),1≤i≤M−1(1.20)u0k=α(tk),uMk=β(tk),0≤k≤N(1.21)
\left\{
\begin{aligned}
&\frac{u_i^{k+1}-u_i^k}{\tau} - a \sdot\frac{1}{2}
\left[
\frac{u_{i+1}^k - 2u_i^k + u_{i-1}^k}{h^2} + \frac{u_{i+1}^{k+1} - 2u_i^{k+1} + u_{i-1}^{k+1}}{h^2}
\right] = f(x_i, t_{k+\frac{\tau}{2}}),\ 1 \le i \le M-1,0 \le k \le N-1 &&(1.19) \\
&u_i^0 = \varphi(x_i), 1 \le i \le M-1 &&(1.20) \\
&u_0^k = \alpha(t_k), u_M^k = \beta(t_k), 0 \le k \le N &&(1.21)
\end{aligned}
\right.
⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧τuik+1−uik−a⋅21[h2ui+1k−2uik+ui−1k+h2ui+1k+1−2uik+1+ui−1k+1]=f(xi,tk+2τ), 1≤i≤M−1,0≤k≤N−1ui0=φ(xi),1≤i≤M−1u0k=α(tk),uMk=β(tk),0≤k≤N(1.19)(1.20)(1.21)
其截断误差如下
Rik(3)=τ224∂3u∂t3(xi,ηik)−ah224[∂4u∂x4(ξik,tk)+∂4u∂x4(ξ~ik,tk+1)]−aτ28∂4u∂x2∂t2(xi,η‾ik),ηik∈(tk,tk+1), η‾ik∈(tk,tk+1), ξik∈(xi−1,xi+1), ξ~ik∈(xi−1,xi+1)
R_{ik}^{(3)} = \frac{\tau^2}{24}\frac{\partial^3 u}{\partial t^3}(x_i, \eta_i^k) - \frac{ah^2}{24}\left[ \frac{\partial^4 u}{\partial x^4}(\xi_i^k, t_k) + \frac{\partial^4 u}{\partial x^4}(\tilde{\xi}_i^k, t_{k+1}) \right] - \frac{a\tau^2}{8}\frac{\partial^4 u}{\partial x^2 \partial t^2}(x_i, \overline{\eta}_i^k),\\
\eta_i^k \in (t_k, t_{k+1}),\ \overline{\eta}_i^k \in (t_k, t_{k+1}),\ \xi_i^k \in (x_{i-1}, x_{i+1}),\ \tilde{\xi}_i^k \in (x_{i-1}, x_{i+1})
Rik(3)=24τ2∂t3∂3u(xi,ηik)−24ah2[∂x4∂4u(ξik,tk)+∂x4∂4u(ξ~ik,tk+1)]−8aτ2∂x2∂t2∂4u(xi,ηik),ηik∈(tk,tk+1), ηik∈(tk,tk+1), ξik∈(xi−1,xi+1), ξ~ik∈(xi−1,xi+1)
差分格式的稳定性和收敛性
- 古典显(隐)格式和 Crank-Nicolson 格式为两层格式,Richardson 格式为三层格式
- 记 Ωh=(x0,x1,...,xM)\Omega_h=(x_0, x_1,..., x_M)Ωh=(x0,x1,...,xM),称 ω=(ω0,ω1,...,ωM)\omega = (\omega_0, \omega_1, ..., \omega_M)ω=(ω0,ω1,...,ωM) 为 Ωh\Omega_hΩh 上的网格函数。记 W={ω∣ωW = \{\omega|\omegaW={ω∣ω 为 Ωh\Omega_hΩh 上的网格函数 }\}},称 WWW 为网格函数空间。
- 设 {uik∣0≤i≤M,0≤k≤N}\{u_i^k|0 \le i \le M,0 \le k \le N \}{uik∣0≤i≤M,0≤k≤N} 是差分格式的解,定义 uk=(u0k,u1k,...,uMk)u^k=(u_0^k, u_1^k,...,u_M^k)uk=(u0k,u1k,...,uMk),则 uku^kuk 为 Ωh\Omega_hΩh 上的网格函数
设 ω∈W\omega \in Wω∈W,定义下面的范数:
L2−范数:∥ω∥2=h[12(w0)2+∑i=1M−1(wi)2+12(ωM)2]L∞−范数:∥ω∥∞=max0≤i≤M∣ωi∣L1−范数:∥ω∥1=h[12∣ω∣+∑i=1M−1∣ωi∣+12∣ωM∣]
\begin{aligned}
&L_2-范数:\parallel\omega\parallel_2 = \sqrt{h\left[ \frac{1}{2}(w_0)^2 + \sum_{i=1}^{M-1}(w_i)^2 + \frac{1}{2}(\omega_M)^2 \right]} \\
&L_{\infty}-范数:\parallel\omega\parallel_{\infty} = \max_{0 \le i \le M}|\omega_i| \\
&L_1-范数:\parallel\omega\parallel_1 = h\left[ \frac{1}{2}|\omega| + \sum_{i=1}^{M-1}|\omega_i| + \frac{1}{2}|\omega_M| \right]
\end{aligned}
L2−范数:∥ω∥2=h[21(w0)2+i=1∑M−1(wi)2+21(ωM)2]L∞−范数:∥ω∥∞=0≤i≤Mmax∣ωi∣L1−范数:∥ω∥1=h[21∣ω∣+i=1∑M−1∣ωi∣+21∣ωM∣]
注:hhh 为空间步长
稳定性
设 {uik∣0≤i≤M,0≤k≤N}\{u_i^k|0 \le i \le M, 0 \le k \le N \}{uik∣0≤i≤M,0≤k≤N} 是差分格式的解,{vik∣0≤i≤M,0≤k≤N}\{v_i^k|0 \le i \le M, 0 \le k \le N \}{vik∣0≤i≤M,0≤k≤N} 是由于初始数据有误差而得到的差分格式的近似解,记
ϵik=uik−vik, 0≤i≤M,0≤k≤N
\epsilon_i^k = u_i^k - v_i^k,\ 0 \le i \le M, 0 \le k \le N
ϵik=uik−vik, 0≤i≤M,0≤k≤N
如果存在与步长 hhh,τ\tauτ 无关的常数 CCC,使得
∥ϵk∥≤C∥ϵ0∥, 1≤k≤N(两层格式)∥ϵk∥≤C(∥ϵ0∥+∥ϵ1∥), 1≤k≤N(三层格式)
\begin{aligned}
&\parallel\epsilon^k\parallel \le C\parallel\epsilon^0\parallel,\ 1 \le k \le N(两层格式) \\
&\parallel\epsilon^k\parallel \le C(\parallel\epsilon^0\parallel+\parallel\epsilon^1\parallel),\ 1 \le k \le N(三层格式)
\end{aligned}
∥ϵk∥≤C∥ϵ0∥, 1≤k≤N(两层格式)∥ϵk∥≤C(∥ϵ0∥+∥ϵ1∥), 1≤k≤N(三层格式)
则称该差分格式关于范数 ∥⋅∥\parallel\sdot\parallel∥⋅∥ 是稳定的,否则称不稳定。
- r=aτh2>0r = \frac{a\tau}{h^2} \gt 0r=h2aτ>0
- 当步长比 r≤12r \le \frac{1}{2}r≤21 时,古典显格式关于 L∞L_{\infty}L∞ 范数是稳定的;当 r>12r \gt \frac{1}{2}r>21 时,其关于 L∞L_{\infty}L∞ 范数不稳定
- 对任意步长比 rrr,古典隐格式关于 L∞L_{\infty}L∞ 范数是稳定的
- 对任意步长比 rrr,Crank-Nicolson 格式关于 L2L_2L2 范数是稳定的
- 对任意步长比 rrr,Richardson 格式关于 L∞L_{\infty}L∞ 范数和 L2L_2L2 范数都是不稳定的
收敛性
设 Uik=u(xi,tk)U_i^k = u(x_i, t_k)Uik=u(xi,tk) 是微分方程定解问题的解,uiku_i^kuik 是对应的差分格式的解,记
eik=Uik−uik,0≤i≤M,0≤k≤N
e_i^k = U_i^k - u_i^k, 0 \le i \le M, 0 \le k \le N
eik=Uik−uik,0≤i≤M,0≤k≤N
若
limh→0,τ→0max0≤k≤N∥ek∥=0
\lim_{h \rightarrow 0,\tau \rightarrow0} \max_{0 \le k \le N}\parallel e^k \parallel = 0
h→0,τ→0lim0≤k≤Nmax∥ek∥=0
则称差分格式在范数 ∥⋅∥\parallel\sdot\parallel∥⋅∥ 下是收敛的。如果
max0≤k≤N∥ek∥=O(hp+τq)
\max_{0 \le k \le N}\parallel e^k \parallel = O(h^p+\tau^q)
0≤k≤Nmax∥ek∥=O(hp+τq)
则称差分格式关于空间步长 ppp 阶、关于时间步长 qqq 阶收敛
- 当步长比 r≤12r \le \frac{1}{2}r≤21 时,古典显格式在 L∞L_{\infty}L∞ 范数下关于空间步长 2 阶、关于时间步长 1 阶收敛
- 对任意步长比 rrr,古典隐格式在 L∞L_{\infty}L∞ 范数下关于空间步长 2 阶、关于时间步长 1 阶收敛
- 对任意步长比 rrr,Crank-Nicolson 格式在 L2L_2L2 范数下关于空间步长和时间步长都是 2 阶收敛的
双曲型方程的差分方法
考虑下面的初边值问题
{∂2u∂t2−a2∂2u∂x2=f(x,t), 0<x<l,0<t<Tu(x,0)=φ(x), ut(x,0)=ψ(x), 0<x<lu(0,t)=α(t), u(l,t)=β(t), 0≤t≤T
\left\{
\begin{aligned}
&\frac{\partial^2u}{\partial t^2} - a^2\frac{\partial^2u}{\partial x^2} = f(x, t),&&\ 0 \lt x \lt l, 0 \lt t \lt T \\
&u(x, 0) = \varphi(x), \ u_t(x, 0) = \psi(x), &&\ 0 \lt x \lt l \\
&u(0, t) = \alpha(t),\ u(l, t) = \beta(t),&&\ 0 \le t \le T
\end{aligned}
\right.
⎩⎪⎪⎪⎨⎪⎪⎪⎧∂t2∂2u−a2∂x2∂2u=f(x,t),u(x,0)=φ(x), ut(x,0)=ψ(x),u(0,t)=α(t), u(l,t)=β(t), 0<x<l,0<t<T 0<x<l 0≤t≤T
记 h=lMh = \frac{l}{M}h=Ml,τ=TN\tau = \frac{T}{N}τ=NT,xi=ihx_i = ihxi=ih,tk=kτt_k=k\tautk=kτ
考虑点 (xi,tk)(x_i, t_k)(xi,tk) 处的方程(显格式和隐格式)
∂2u∂t2(xi,tk)−a2∂2u∂x2(xi,tk)=f(xi,tk)
\frac{\partial^2 u}{\partial t^2}(x_i, t_k) - a^2\frac{\partial^2 u}{\partial x^2}(x_i, t_k) = f(x_i, t_k)
∂t2∂2u(xi,tk)−a2∂x2∂2u(xi,tk)=f(xi,tk)
- 利用二阶差商 (0.4)(0.4)(0.4)、初边值条件和泰勒展开可得显格式
- 利用二阶差商 (0.4)(0.4)(0.4)、平均公式 (0.5)(0.5)(0.5) 和初边值条件可得隐格式
显格式
{1τ2(uik+1−2uik+uik−1)−a2h2(ui+1k−2uik+ui−1k)=f(xi,tk), 1≤i≤M−1,1≤k≤N−1(3.2)ui0=φ(xi), ui1=Ψ(xi), 1≤i≤M−1(3.3)u0k=α(tk), uMk=β(tk), 0<k<N(3.4) \left\{ \begin{aligned} &\frac{1}{\tau^2}(u_i^{k+1} - 2u_i^k + u_i^{k-1}) - \frac{a^2}{h^2}(u_{i+1}^k - 2u_i^k + u_{i-1}^k) = f(x_i, t_k), \ 1 \le i \le M-1,1 \le k \le N-1 &&(3.2)\\ &u_i^0 = \varphi(x_i), \ u_i^1 = \Psi(x_i), \ 1 \le i \le M-1 &&(3.3)\\ &u_0^k = \alpha(t_k), \ u_M^k = \beta(t_k), \ 0 \lt k \lt N &&(3.4) \end{aligned} \right. ⎩⎪⎪⎪⎨⎪⎪⎪⎧τ21(uik+1−2uik+uik−1)−h2a2(ui+1k−2uik+ui−1k)=f(xi,tk), 1≤i≤M−1,1≤k≤N−1ui0=φ(xi), ui1=Ψ(xi), 1≤i≤M−1u0k=α(tk), uMk=β(tk), 0<k<N(3.2)(3.3)(3.4)
注:Ψ(xi)=φ(xi)+τψ(xi)+τ22[a2φ′′(xi)+f(xi,0)]\Psi(x_i) = \varphi(x_i)+\tau\psi(x_i)+\frac{\tau^2}{2}[a^2\varphi^{''}(x_i)+f(x_i,0)]Ψ(xi)=φ(xi)+τψ(xi)+2τ2[a2φ′′(xi)+f(xi,0)]
其截断误差如下
Rik=τ212∂4u∂t4(xi,ηik)−a2h212∂4u∂x4(ξik,tk), ηik∈(tk−1,tk+1), ξik∈(xi−1,xi+1)
R_{ik} = \frac{\tau^2}{12}\frac{\partial^4 u}{\partial t^4}(x_i, \eta_i^k) - \frac{a^2h^2}{12}\frac{\partial^4 u}{\partial x^4}(\xi_i^k, t_k),\ \eta_i^k \in (t_{k-1}, t_{k+1}),\ \xi_i^k \in (x_{i-1}, x_{i+1})
Rik=12τ2∂t4∂4u(xi,ηik)−12a2h2∂x4∂4u(ξik,tk), ηik∈(tk−1,tk+1), ξik∈(xi−1,xi+1)
记 s=aτhs = \frac{a\tau}{h}s=haτ,sss 称为步长比,则可将上面差分格式改写为
uik+1=s2(ui+1k+ui−1k)+2(1−s2)uik−uik−1+τ2f(xi,tk), 1≤i≤M−1,1≤k≤N−1
u_i^{k+1} = s^2(u_{i+1}^k + u_{i-1}^k) + 2(1-s^2)u_i^k - u_i^{k-1} + \tau^2f(x_i, t_k),\ 1 \le i \le M-1,1 \le k \le N-1 \\
uik+1=s2(ui+1k+ui−1k)+2(1−s2)uik−uik−1+τ2f(xi,tk), 1≤i≤M−1,1≤k≤N−1
- 当步长比 s≤1s \le 1s≤1 时,显式差分格式 (3.2)−(3.4)(3.2)-(3.4)(3.2)−(3.4) 在 L2L_2L2 范数下是稳定的;当步长比 s>1s \gt 1s>1 时,在 L2L_2L2 范数下是不稳定的
- 当步长比 s≤1s \le 1s≤1 时,显式差分格式 (3.2)−(3.4)(3.2)-(3.4)(3.2)−(3.4) 在 L2L_2L2 范数下关于空间步长和时间步长都是二阶收敛的
隐格式
{1τ2(uik+1−2uik+uik−1)−a22h2(ui+1k+1−2uik+1+ui−1k+1+ui+1k−1−2uik−1+ui−1k−1)=f(xi,tk), 1≤i≤M−1,1≤k≤N−1(3.6)ui0=φ(xi), ui1=Ψ(x1), 1≤i≤M−1(3.7)u0k=α(tk),uMk=β(tk), 0≤k≤N(3.8) \left\{ \begin{aligned} \frac{1}{\tau^2}&(u_i^{k+1} - 2u_i^k + u_i^{k-1}) - \frac{a^2}{2h^2}(u_{i+1}^{k+1} - 2u_i^{k+1} +u_{i-1}^{k+1} + u_{i+1}^{k-1} - 2u_i^{k-1} + u_{i-1}^{k-1}) \\ &= f(x_i, t_k), \ 1 \le i \le M-1,1 \le k \le N-1 && (3.6)\\ u_i^0 &= \varphi(x_i), \ u_i^1 = \Psi(x_1), \ 1 \le i \le M-1 && (3.7) \\ u_0^k &= \alpha(t_k), u_M^k = \beta(t_k), \ 0 \le k \le N && (3.8) \end{aligned} \right. ⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧τ21ui0u0k(uik+1−2uik+uik−1)−2h2a2(ui+1k+1−2uik+1+ui−1k+1+ui+1k−1−2uik−1+ui−1k−1)=f(xi,tk), 1≤i≤M−1,1≤k≤N−1=φ(xi), ui1=Ψ(x1), 1≤i≤M−1=α(tk),uMk=β(tk), 0≤k≤N(3.6)(3.7)(3.8)
- 对任意步长比 sss,隐式差分格式 (3.6)−(3.8)(3.6)-(3.8)(3.6)−(3.8) 在 L2L_2L2 范数下是稳定的
- 对任意步长比 sss,上述隐式差分格式在 L2L_2L2 范数下关于空间步长和时间步长都是二阶收敛的
其截断误差如下
Rik=−a2τ22∂4u∂x2∂t2(xi,η‾ik)+τ212∂4u∂t4(xi,ηik)−a2h224[∂4u∂x4(ξ‾ik+1,tk+1)+∂4u∂x4(ξ‾ik−1,tk−1)]ηik∈(tk−1,tk+1), η‾ik∈(tk−1,tk+1), ξ‾i∈(xi−1,xi+1)
R_{ik} = -\frac{a^2\tau^2}{2}\frac{\partial^4 u}{\partial x^2 \partial t^2}(x_i, \overline{\eta}_i^k) + \frac{\tau^2}{12}\frac{\partial^4 u}{\partial t^4}(x_i, \eta_i^k) - \frac{a^2 h^2}{24}\left[ \frac{\partial^4 u}{\partial x^4}(\overline{\xi}_i^{k+1}, t_{k+1}) + \frac{\partial^4 u}{\partial x^4}(\overline{\xi}_i^{k-1}, t_{k-1}) \right] \\
\eta_i^k \in (t_{k-1}, t_{k+1}),\ \overline{\eta}_i^k \in (t_{k-1}, t_{k+1}), \ \overline{\xi}_i \in (x_{i-1}, x_{i+1})
Rik=−2a2τ2∂x2∂t2∂4u(xi,ηik)+12τ2∂t4∂4u(xi,ηik)−24a2h2[∂x4∂4u(ξik+1,tk+1)+∂x4∂4u(ξik−1,tk−1)]ηik∈(tk−1,tk+1), ηik∈(tk−1,tk+1), ξi∈(xi−1,xi+1)
椭圆型方程的差分方法
考虑二维 Poisson 方程 Dirichlet 边值问题
{−(∂2u∂x2+∂2u∂y2)=f(x,y)(x,y)∈Ω(4.1)u=φ(x,y)(x,y)∈Γ(4.2)
\left\{
\begin{aligned}
&-(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}) = f(x, y) && (x, y) \in \Omega && (4.1) \\
&u = \varphi(x, y) && (x, y) \in \Gamma && (4.2)
\end{aligned}
\right.
⎩⎪⎨⎪⎧−(∂x2∂2u+∂y2∂2u)=f(x,y)u=φ(x,y)(x,y)∈Ω(x,y)∈Γ(4.1)(4.2)
其中 Ω\OmegaΩ 为矩形区域 Ω={(x,y)∣a<x<b,c<y<d}\Omega = \{(x, y)|a \lt x \lt b, c \lt y \lt d \}Ω={(x,y)∣a<x<b,c<y<d},Γ\GammaΓ 是 Ω\OmegaΩ 的边界
差分格式的建立
将区间 [a,b][a, b][a,b] 作 mmm 等分,将 [c,d][c, d][c,d] 作 nnn 等分,记 h1=b−amh_1 = \frac{b-a}{m}h1=mb−a,h2=d−cnh_2 = \frac{d-c}{n}h2=nd−c,xi=a+ih1(0≤i≤m)x_i = a + ih_1(0 \le i \le m)xi=a+ih1(0≤i≤m),yj=c+jh2(0≤j≤n)y_j = c + jh_2(0 \le j \le n)yj=c+jh2(0≤j≤n)。称 h1h_1h1 和 h2h_2h2 为 xxx 方向和 yyy 方向的步长,用平行线
x=xi0≤i≤my=yj0≤j≤n
\begin{aligned}
&x = x_i && 0 \le i \le m \\
&y = y_j && 0 \le j \le n
\end{aligned}
x=xiy=yj0≤i≤m0≤j≤n
将 Ω\OmegaΩ 剖分为 mnmnmn 个小矩形,交点 (xi,yj)(x_i, y_j)(xi,yj) 称为网格节点,记
Ωh={(xi,yj)∣0≤i≤m,0≤j≤n}Ωh∘={(xi,yj)∣1≤i≤m−1,1≤j≤n−1}
\begin{aligned}
&\Omega_h = \{ (x_i, y_j)| 0 \le i \le m, 0 \le j \le n \} \\
&\Omega_h^{\circ} = \{ (x_i, y_j)| 1 \le i \le m-1, 1 \le j \le n-1 \} \\
\end{aligned}
Ωh={(xi,yj)∣0≤i≤m,0≤j≤n}Ωh∘={(xi,yj)∣1≤i≤m−1,1≤j≤n−1}
称 Ωh∘\Omega_h^{\circ}Ωh∘ 为内节点,称 Γh=Ωh\Gamma_h = \Omega_hΓh=Ωh\ Ωh∘\Omega_h^{\circ}Ωh∘ 为边界节点。记
ω≡{(i,j)∣(xi,yj)∈Ωh∘},γ≡{(i,j)∣(xi,yj)∈Γh}
\omega \equiv \{ (i, j)|(x_i, y_j) \in \Omega_h^{\circ} \},\quad \gamma \equiv \{ (i, j)|(x_i, y_j) \in \Gamma_h \}
ω≡{(i,j)∣(xi,yj)∈Ωh∘},γ≡{(i,j)∣(xi,yj)∈Γh}
在节点 (xi,yj)(x_i, y_j)(xi,yj) 处考虑边值问题
{−[∂2u∂x2(xi,yj)+∂2u∂y2(xi,yj)]=f(xi,yj)(i,j)∈ω(4.3)u(xi,yj)=φ(xi,yj)(i,j)∈γ(4.4)
\left\{
\begin{aligned}
&-\left[\frac{\partial^2 u}{\partial x^2}(x_i, y_j) + \frac{\partial^2 u}{\partial y^2}(x_i, y_j)\right] = f(x_i, y_j) && (i, j) \in \omega && (4.3) \\
&u(x_i, y_j) = \varphi(x_i, y_j) && (i, j) \in \gamma && (4.4)
\end{aligned}
\right.
⎩⎪⎨⎪⎧−[∂x2∂2u(xi,yj)+∂y2∂2u(xi,yj)]=f(xi,yj)u(xi,yj)=φ(xi,yj)(i,j)∈ω(i,j)∈γ(4.3)(4.4)
五点差分格式
{−(δx2u+δy2u)=f(xi,yj)(i,j)∈ω(4.8)uij=φ(xi,yj)(i,j)∈γ(4.9)δx2u=1h12[ui−1,j−2uij+ui+1,j]δy2u=1h22[ui,j−1−2uij+ui,j+1] \left\{ \begin{aligned} &-(\delta_x^2 u + \delta_y^2 u) = f(x_i, y_j) && (i, j) \in \omega && (4.8) \\ &u_{ij} = \varphi(x_i, y_j) && (i, j) \in \gamma && (4.9) \end{aligned} \right. \\ \\ \delta_x^2u = \frac{1}{h_1^2} \left[ u_{i-1,j} - 2u_{ij} + u_{i+1, j} \right] \\ \delta_y^2u = \frac{1}{h_2^2} \left[ u_{i,j-1} - 2u_{ij} + u_{i, j+1} \right] \\ {−(δx2u+δy2u)=f(xi,yj)uij=φ(xi,yj)(i,j)∈ω(i,j)∈γ(4.8)(4.9)δx2u=h121[ui−1,j−2uij+ui+1,j]δy2u=h221[ui,j−1−2uij+ui,j+1]
注:Uij=u(xi,yj)U_{ij} = u(x_i, y_j)Uij=u(xi,yj),而 uiju_{ij}uij 为 UijU_{ij}Uij 的近似
其截断误差如下
Rij=−1h12[u(xi−1,yj)−2u(xi,yj)+u(xi+1,yj)]−1h22[u(xi,yj−1)−2u(xi,yj)+u(xi,yj+1)]−f(xi,yj)
R_{ij} = -\frac{1}{h_1^2}\left[ u(x_{i-1}, y_j) - 2u(x_i, y_j) + u(x_{i+1}, y_j) \right] - \frac{1}{h_2^2}\left[ u(x_i, y_{j-1}) - 2u(x_i, y_j) + u(x_i, y_{j+1}) \right] - f(x_i, y_j)
Rij=−h121[u(xi−1,yj)−2u(xi,yj)+u(xi+1,yj)]−h221[u(xi,yj−1)−2u(xi,yj)+u(xi,yj+1)]−f(xi,yj)
- 差分格式 (4.8)−(4.9)(4.8)-(4.9)(4.8)−(4.9) 存在唯一解
- 设 {u(x,y)∣a≤x≤b,c≤y≤d}\{ u(x, y)|a \le x \le b, c \le y \le d \}{u(x,y)∣a≤x≤b,c≤y≤d} 是定解问题 (4.1)−(4.2)(4.1)-(4.2)(4.1)−(4.2) 的解,{uij∣0≤i≤m,0≤j≤n}\{ u_{ij}| 0 \le i \le m, 0 \le j \le n \}{uij∣0≤i≤m,0≤j≤n} 为差分格式 (4.8)−(4.9)(4.8)-(4.9)(4.8)−(4.9) 的解,则有
max(i,j)∈ω∣u(xi,yj)−uij∣≤M448[(b−a2)2+(d−c2)2](h12+h22)M4=max{max(x,y)∈Ω‾∣∂4u(x,y)∂x4∣,max(x,y)∈Ω‾∣∂4u(x,y)∂y4∣} \begin{aligned} &\max_{(i, j) \in \omega} \big| u(x_i, y_j) - u_{ij} \big| \le \frac{M_4}{48} \left[ (\frac{b-a}{2})^2 + (\frac{d-c}{2})^2 \right](h_1^2 + h_2^2) \\\\ &M_4 = \max\left\{ \max_{(x, y) \in \overline{\Omega}}\bigg|\frac{\partial^4u(x, y)}{\partial x^4}\bigg|, \max_{(x, y) \in \overline{\Omega}}\bigg|\frac{\partial^4u(x, y)}{\partial y^4}\bigg| \right\} \end{aligned} (i,j)∈ωmax∣∣u(xi,yj)−uij∣∣≤48M4[(2b−a)2+(2d−c)2](h12+h22)M4=max{(x,y)∈Ωmax∣∣∣∣∂x4∂4u(x,y)∣∣∣∣,(x,y)∈Ωmax∣∣∣∣∂y4∂4u(x,y)∣∣∣∣}
- 当步长 h1h_1h1 和 h2h_2h2 同时缩小到原来的 12\frac{1}{2}21 时,最大误差约缩小到原来的 14\frac{1}{4}41
4815

被折叠的 条评论
为什么被折叠?



