PDE抛物型方程數值解法總結與例題分析

抛物型方程

例題及解答

例題：

構造抛物型方程

{ ∂ u ∂ t = ∂ ∂ x ( x ∂ u ∂ x ) , 0.5 < x < 1 , 0 < t ⩽ T , u ( x ,   0 ) = φ ( x ) , 0.5 ⩽ x ⩽ 1 , u ( 0.5 ,   t ) = 0 ,   ∂ u ∂ x ( 1 , t ) = − 1 3 u ( 1 ,   t ) ,   0 ⩽ t ⩽ T , \left\{\begin{aligned} &\frac{\partial u}{\partial t}=\frac{\partial}{\partial x}\left(x\frac{\partial u}{\partial x}\right),\quad0.5<x<1,0<t\leqslant T,\\ &u(x, \ 0)=\varphi(x),\quad0.5\leqslant x \leqslant 1,\\ &u(0.5,\ t)=0,\ \frac{\partial u}{\partial x}(1,t)=-\frac13 u(1,\ t), \ 0\leqslant t \leqslant T, \end{aligned}\right. ⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧​​∂t∂u​=∂x∂​(x∂x∂u​),0.5<x<1,0<t⩽T,u(x, 0)=φ(x),0.5⩽x⩽1,u(0.5, t)=0, ∂x∂u​(1,t)=−31​u(1, t), 0⩽t⩽T,​

的顯格式，取 h = 0.1 ,   r = 0.1 h=0.1,\ r=0.1 h=0.1, r=0.1。

解答：

由于本題空間項同時含有一階項和二階項，導緻截斷誤差的階數不比對（應該同為二階），利用有限差分格式不易求解，故下面用積分插值進行求解。

用積分插值法建立差分格式如下：

令 x = j h ,   t = k τ ,   j = 0 , 1 , ⋯   , 1 / h , k = 0 , 1 , ⋯   , T / τ x=jh,\ t=k\tau,\ j=0,1,\cdots,1/h,\quad k=0,1,\cdots,T/\tau x=jh, t=kτ, j=0,1,⋯,1/h,k=0,1,⋯,T/τ

設 ( j , k ) (j,k) (j,k)為内點，取區間 x j − h 2 ⩽ x ⩽ x j + h 2 x_j-\frac h2\leqslant x\leqslant x_j+\frac h2 xj​−2h​⩽x⩽xj​+2h​，在此區間對 x x x積分得：

∫ x j − 1 2 x j + 1 2 ∂ u ∂ t d x = ∫ x j − 1 2 x j + 1 2 ∂ ∂ x ( x ∂ u ∂ x ) d x , \int_{x_{j-\frac 12}}^{x_{j+\frac 12}}\frac{\partial u}{\partial t}\mathrm{d}x= \int_{x_{j-\frac 12}}^{x_{j+\frac 12}}\frac{\partial }{\partial x}\left( x\frac{\partial u}{\partial x}\right) \mathrm{d}x, ∫xj−21​​xj+21​​​∂t∂u​dx=∫xj−21​​xj+21​​​∂x∂​(x∂x∂u​)dx,

利用數值積分的中矩形公式得：

∫ x j − 1 2 x j + 1 2 ∂ u ∂ t d x = [ ∂ u ∂ t ] j h + O ( h 3 ) , \int_{x_{j-\frac 12}}^{x_{j+\frac 12}}\frac{\partial u}{\partial t}\mathrm{d}x=\left[ \frac{\partial u}{\partial t}\right]_j h +\mathcal{O}(h^3) , ∫xj−21​​xj+21​​​∂t∂u​dx=[∂t∂u​]j​h+O(h3),

右邊同理：

∫ x j − 1 2 x j + 1 2 ∂ ∂ x ( x ∂ u ∂ x ) d x = x j + 1 2 [ ∂ u ∂ x ] j + 1 2 − x j − 1 2 [ ∂ u ∂ x ] j − 1 2 , (1) \int_{x_{j-\frac 12}}^{x_{j+\frac 12}}\frac{\partial }{\partial x}\left( x\frac{\partial u}{\partial x}\right) \mathrm{d}x = x_{j+\frac12}\left[ \frac{\partial u}{\partial x}\right]_{j+\frac12} -x_{j-\frac12}\left[ \frac{\partial u}{\partial x}\right]_{j-\frac12}, \tag{1} ∫xj−21​​xj+21​​​∂x∂​(x∂x∂u​)dx=xj+21​​[∂x∂u​]j+21​​−xj−21​​[∂x∂u​]j−21​​,(1)

由最簡顯格式公式，可得

[ ∂ u ∂ t ] j = [ u ] j k + 1 − [ u ] j k τ + O ( τ ) . \left[ \frac{\partial u}{\partial t}\right]_j=\frac{[u]_j^{k+1}-[u]_j^k}{\tau}+\mathcal{O}(\tau). [∂t∂u​]j​=τ[u]jk+1​−[u]jk​​+O(τ).

對式 ( 1 ) (1) (1) 右端用一階中心差商代替一階偏導數得到

[ ∂ u ∂ x ] j + 1 2 = [ u ] j + 1 k − [ u ] j k h + O ( h 2 ) , \left[ \frac{\partial u}{\partial x}\right]_{j+\frac12}=\frac{[u]_{j+1}^{k}-[u]_j^k}{h}+\mathcal{O}(h^2), [∂x∂u​]j+21​​=h[u]j+1k​−[u]jk​​+O(h2),

[ ∂ u ∂ x ] j − 1 2 = [ u ] j k − [ u ] j − 1 k h + O ( h 2 ) , \left[ \frac{\partial u}{\partial x}\right]_{j-\frac12}=\frac{[u]_{j}^{k}-[u]_{j-1}^k}{h}+\mathcal{O}(h^2), [∂x∂u​]j−21​​=h[u]jk​−[u]j−1k​​+O(h2),

聯立上面各式，得到：

[ u ] j k + 1 − [ u ] j k τ = 1 h 2 [ x j + 1 2 ( [ u ] j + 1 k − [ u ] j k ) − x j − 1 2 ( [ u ] j + 1 k − [ u ] j − 1 k ) ] + O ( τ + h 2 ) \frac{[u]_j^{k+1}-[u]_j^k}{\tau}=\frac{1}{h^2}\left[ x_{j+\frac12}\left( [u]_{j+1}^{k}-[u]_j^k\right) -x_{j-\frac12}\left( [u]_{j+1}^{k}-[u]_{j-1}^k\right) \right] +\mathcal{O}(\tau+h^2) τ[u]jk+1​−[u]jk​​=h21​[xj+21​​([u]j+1k​−[u]jk​)−xj−21​​([u]j+1k​−[u]j−1k​)]+O(τ+h2)

用 u j k u_j^k ujk​代替 [ u ] j k [u]_j^k [u]jk​，去掉誤差項，則有

u j k + 1 − u j k τ = 1 h 2 [ x j + 1 2 ( u j + 1 k − u j k ) − x j − 1 2 ( u j k − u j − 1 k ) ] \frac{u_j^{k+1}-u_j^k}{\tau}=\frac{1}{h^2}\left[ x_{j+\frac12}\left( u_{j+1}^{k}-u_j^k\right) -x_{j-\frac12}\left( u_{j}^{k}-u_{j-1}^k\right) \right] τujk+1​−ujk​​=h21​[xj+21​​(uj+1k​−ujk​)−xj−21​​(ujk​−uj−1k​)]

取網比為 r = τ h 2 r=\frac \tau {h^2} r=h2τ​，則上式等價于

u j k + 1 = r x j − 1 2 u j − 1 k + ( 1 − r x j + 1 2 − r x j − 1 2 ) u j k + r x j + 1 2 u j + 1 k . u_j^{k+1}=rx_{j-\frac12}u_{j-1}^k+(1-rx_{j+\frac12}-rx_{j-\frac12})u_j^k+rx_{j+\frac12}u_{j+1}^k. ujk+1​=rxj−21​​uj−1k​+(1−rxj+21​​−rxj−21​​)ujk​+rxj+21​​uj+1k​.

代入步長得：

u j k + 1 = 0.1 x j − 1 2 u j − 1 k + ( 1 − 0.1 x j + 1 2 − 0.1 x j − 1 2 ) u j k + 0.1 x j + 1 2 u j + 1 k . u_j^{k+1}=0.1x_{j-\frac12}u_{j-1}^k+(1-0.1x_{j+\frac12}-0.1x_{j-\frac12})u_j^k+0.1x_{j+\frac12}u_{j+1}^k. ujk+1​=0.1xj−21​​uj−1k​+(1−0.1xj+21​​−0.1xj−21​​)ujk​+0.1xj+21​​uj+1k​.

而由于 0.5 < x < 1 0.5<x<1 0.5<x<1，是以 x j > 0 ,   x j − 1 2 > 0 ,   x j + 1 2 > 0 ,   j = 0 , 1 , ⋯   , 1 / h . x_j>0,\ x_{j-\frac12}>0,\ x_{j+\frac12}>0,\ j=0,1,\cdots,1/h. xj​>0, xj−21​​>0, xj+21​​>0, j=0,1,⋯,1/h.

對第一類邊界條件 u ( x , 0 ) = φ ( x ) ,    u ( 0.5 , t ) = 0 u(x,0)=\varphi(x),\ \ u(0.5,t)=0 u(x,0)=φ(x),  u(0.5,t)=0，有：

u j 0 = φ ( x j ) ,    u 5 k = 0 , u_j^0=\varphi(x_j),\ \ u_{5}^k=0, uj0​=φ(xj​),  u5k​=0,

對第二類邊界條件 ∂ ∂ x u ( 1 , t ) = − 1 3 u ( 1 ,   t ) \frac{\partial }{\partial x}u(1,t)=-\frac13 u(1,\ t) ∂x∂​u(1,t)=−31​u(1, t)，可得到：

[ ∂ u ∂ x ] j k = u j + 1 k − u j − 1 k 2 h , \left[ \frac{\partial u}{\partial x}\right]_j^k=\frac{u_{j+1}^{k}-u_{j-1}^k}{2h}, [∂x∂u​]jk​=2huj+1k​−uj−1k​​,

是以

[ ∂ u ∂ x ] 10 k = u 11 k − u 9 k 0.2 = − 1 3 u 10 k , \left[ \frac{\partial u}{\partial x}\right]_{10}^k=\frac{u_{11}^{k}-u_{9}^k}{0.2}=-\frac13 u_{10}^k, [∂x∂u​]10k​=0.2u11k​−u9k​​=−31​u10k​,

即： u 10 k = 15 ( u 9 k − u 11 k ) u_{10}^k=15(u_9^k-u_{11}^k) u10k​=15(u9k​−u11k​).

小結

本文總結了一道特殊形式的抛物型方程的數值解格式，其需要掌握的知識主要是有限差分格式，即不同的差商表示對偏導項的離散，其難點在于對數值分析方法的掌握，同時還要注意方程的自變量範圍對離散格式的影響以及初值邊值條件的處理。

參考

[1] 常微分數值方法部分課後習題答案