数值分析复习(七)——偏微分方程数值解法

七、偏微分方程数值解法

偏微分方程一般可分为抛物型方程、双曲型方程和椭圆型方程三种类型,其一般形式如下
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=B24AC=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(x0h)]+2hg(ξ2), ξ2(x0h,x0)g(x0)=h1[g(x0+2h)g(x02h)]24h2g(ξ3), ξ3(x02h,x0+2h)g(x0)=h21[g(x0+h)2g(x0)+g(x0h)]12h2g(4)(ξ4), ξ4(x0h,x0+h)g(x0)=21[g(x0h)+g(x0+h)]2h2g(ξ5), ξ5(x0h,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. tuax22u=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)0xl,0tT}h=lMh=\frac{l}{M}h=Mlτ=TN\tau = \frac{T}{N}τ=NTxi=ih(0≤i≤M)x_i=ih(0 \le i \le M)xi=ih(0iM)tk=kτ(0≤k≤N)t_k = k\tau(0 \le k \le N)tk=kτ(0kN),用两簇平行直线
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)0iM,0kN}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)1iM1,1kN}Γ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) tu(xi,tk)ax22u(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+1uik)h2a(ui+1k2uik+ui1k)=f(xi,tk), 1iM1,0kN1ui0=φ(xi), 1iM1u0k=α(tk), uMk=β(tk), 0kN(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=(12r)uik+r(ui1k+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+1uM2k+1uM1k+1=12rrr12rrr12rrr12ru1ku2kuM2kuM1k+τf(x1,tk)+rα(tk)τf(x2,tk)τf(xM2,tk)τf(xM1,tk)+rβ(tk)k=0,1,...,N1

其截断误差如下
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τt22u(xi,ηik)12ah2x44u(ξik,tk), ηik(tk,tk+1), ξik(xi1,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(uikuik1)h2a(ui+1k2uik+ui1k)=f(xi,tk), 1iM1,1kNui0=φ(xi), 1iM1u0k=α(tk), uMk=β(tk), 0kN(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+2rrr1+2rrr1+2rrr1+2ru1ku2kuM2kuM1k=u1k1u2k1uM2k1uM1k1+τf(x1,tk)+rα(tk)τf(x2,tk)τf(xM2,tk)τf(xM1,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τt22u(xi,ηik)12ah2x44u(ξik,tk), ηik(tk1,tk), ξik(xi1,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}}) tu(xi,tk+2τ)ax22u(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+1uika21[h2ui+1k2uik+ui1k+h2ui+1k+12uik+1+ui1k+1]=f(xi,tk+2τ), 1iM1,0kN1ui0=φ(xi),1iM1u0k=α(tk),uMk=β(tk),0kN(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τ2t33u(xi,ηik)24ah2[x44u(ξik,tk)+x44u(ξ~ik,tk+1)]8aτ2x2t24u(xi,ηik),ηik(tk,tk+1), ηik(tk,tk+1), ξik(xi1,xi+1), ξ~ik(xi1,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 \}{uik0iM,0kN} 是差分格式的解,定义 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∞−范数:∥ω∥∞=max⁡0≤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=1M1(wi)2+21(ωM)2]Lω=0iMmaxωiL1ω1=h[21ω+i=1M1ωi+21ωM]
hhh 为空间步长

稳定性

{uik∣0≤i≤M,0≤k≤N}\{u_i^k|0 \le i \le M, 0 \le k \le N \}{uik0iM,0kN} 是差分格式的解,{vik∣0≤i≤M,0≤k≤N}\{v_i^k|0 \le i \le M, 0 \le k \le N \}{vik0iM,0kN} 是由于初始数据有误差而得到的差分格式的近似解,记
ϵ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=uikvik, 0iM,0kN
如果存在与步长 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} ϵkCϵ0, 1kNϵkC(ϵ0+ϵ1), 1kN
则称该差分格式关于范数 ∥⋅∥\parallel\sdot\parallel 是稳定的,否则称不稳定。

  • r=aτh2>0r = \frac{a\tau}{h^2} \gt 0r=h2aτ>0
  • 当步长比 r≤12r \le \frac{1}{2}r21 时,古典显格式关于 L∞L_{\infty}L 范数是稳定的;当 r>12r \gt \frac{1}{2}r>21 时,其关于 L∞L_{\infty}L 范数不稳定
  • 对任意步长比 rrr古典隐格式关于 L∞L_{\infty}L 范数是稳定的
  • 对任意步长比 rrrCrank-Nicolson 格式关于 L2L_2L2 范数是稳定的
  • 对任意步长比 rrrRichardson 格式关于 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=Uikuik,0iM,0kN

lim⁡h→0,τ→0max⁡0≤k≤N∥ek∥=0 \lim_{h \rightarrow 0,\tau \rightarrow0} \max_{0 \le k \le N}\parallel e^k \parallel = 0 h0,τ0lim0kNmaxek=0
则称差分格式在范数 ∥⋅∥\parallel\sdot\parallel 下是收敛的。如果
max⁡0≤k≤N∥ek∥=O(hp+τq) \max_{0 \le k \le N}\parallel e^k \parallel = O(h^p+\tau^q) 0kNmaxek=O(hp+τq)
则称差分格式关于空间步长 ppp 阶、关于时间步长 qqq 阶收敛

  • 当步长比 r≤12r \le \frac{1}{2}r21 时,古典显格式在 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. t22ua2x22u=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 0tT
h=lMh = \frac{l}{M}h=Mlτ=TN\tau = \frac{T}{N}τ=NTxi=ihx_i = ihxi=ihtk=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) t22u(xi,tk)a2x22u(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+12uik+uik1)h2a2(ui+1k2uik+ui1k)=f(xi,tk), 1iM1,1kN1ui0=φ(xi), ui1=Ψ(xi), 1iM1u0k=α(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τ2t44u(xi,ηik)12a2h2x44u(ξik,tk), ηik(tk1,tk+1), ξik(xi1,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+ui1k)+2(1s2)uikuik1+τ2f(xi,tk), 1iM1,1kN1

  • 当步长比 s≤1s \le 1s1 时,显式差分格式 (3.2)−(3.4)(3.2)-(3.4)(3.2)(3.4)L2L_2L2 范数下是稳定的;当步长比 s>1s \gt 1s>1 时,在 L2L_2L2 范数下是不稳定
  • 当步长比 s≤1s \le 1s1 时,显式差分格式 (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+12uik+uik1)2h2a2(ui+1k+12uik+1+ui1k+1+ui+1k12uik1+ui1k1)=f(xi,tk), 1iM1,1kN1=φ(xi), ui1=Ψ(x1), 1iM1=α(tk),uMk=β(tk), 0kN(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τ2x2t24u(xi,ηik)+12τ2t44u(xi,ηik)24a2h2[x44u(ξik+1,tk+1)+x44u(ξik1,tk1)]ηik(tk1,tk+1), ηik(tk1,tk+1), ξi(xi1,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. (x22u+y22u)=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=mbah2=d−cnh_2 = \frac{d-c}{n}h2=ndcxi=a+ih1(0≤i≤m)x_i = a + ih_1(0 \le i \le m)xi=a+ih1(0im)yj=c+jh2(0≤j≤n)y_j = c + jh_2(0 \le j \le n)yj=c+jh2(0jn)。称 h1h_1h1h2h_2h2xxx 方向和 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=yj0im0jn
Ω\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)0im,0jn}Ωh={(xi,yj)1im1,1jn1}
Ω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. [x22u(xi,yj)+y22u(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[ui1,j2uij+ui+1,j]δy2u=h221[ui,j12uij+ui,j+1]

Uij=u(xi,yj)U_{ij} = u(x_i, y_j)Uij=u(xi,yj),而 uiju_{ij}uijUijU_{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(xi1,yj)2u(xi,yj)+u(xi+1,yj)]h221[u(xi,yj1)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)axb,cyd} 是定解问题 (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 \}{uij0im,0jn} 为差分格式 (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)ωmaxu(xi,yj)uij48M4[(2ba)2+(2dc)2](h12+h22)M4=max{(x,y)Ωmaxx44u(x,y),(x,y)Ωmaxy44u(x,y)}

  • 当步长 h1h_1h1h2h_2h2 同时缩小到原来的 12\frac{1}{2}21 时,最大误差约缩小到原来的 14\frac{1}{4}41
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值