logo

ディリクレ境界条件が与えられた熱方程式の初期値問題に対する数値解析的解法 📂数値解析

ディリクレ境界条件が与えられた熱方程式の初期値問題に対する数値解析的解法

##例

{ut=γuxxu(t,0)=u(t,l)=0u(0,x)=f(x) \begin{cases} u_{t} = \gamma u_{xx} \\ u(t,0) = u(t,l) = 0 \\ u(0,x) = f(x) \end{cases}

与えられた問題には代数的解があるほど簡単で、しかし、微分方程式を解くための数値解析法をなぜ学ぶのかを明快に示す例にもなる。これは単にy′=f(x,y)y ' = f(x,y)形式の微分方程式を解くことが、偏微分方程式の解法にも繋がるからだ。

解法

代数的解法とは異なり、γ\gammaの値は重要ではないので、便宜上γ=1\gamma = 1としよう。

xj=δj x_{j} = \delta j 直線空間line spacexxのステップサイズをδ\deltaとすると、j=0,1,⋯ ,mj= 0, 1, \cdots , mに対して上記のように表せる。

xxに対する二階微分の近似として Uxx≈U(t,xj+1)−2U(t,xj)+U(t,xj−1)δ2 U_{xx} \approx {{ U(t, x_{j+1}) - 2 U(t, x_{j}) + U(t, x_{j-1}) } \over { \delta^2 }} これでttに関する方程式 ut(t,xj)=u(t,xj+1)−2u(t,xj)+u(t,xj−1)δ2 u_{t} (t, x_{j}) = {{ u(t, x_{j+1}) - 2 u(t, x_{j}) + u(t, x_{j-1}) } \over { \delta^2 }} をm−1m-1個得ることができる。これを行列で表すと [ut(t,x1)ut(t,x2)ut(t,x3)⋮ut(t,xm−2)ut(t,xm−1)]=1δ2[−210⋯001−21⋯0001−2⋯00⋮⋮⋮⋱⋮⋮000⋯−21000⋯1−2][u(t,x1)u(t,x2)u(t,x3)⋮u(t,xm−2)u(t,xm−1)] \begin{bmatrix} u_{t} (t, x_{1}) \\ u_{t} (t, x_{2}) \\ u_{t} (t, x_{3}) \\ \vdots \\ u_{t} (t, x_{m-2}) \\ u_{t} (t, x_{m-1}) \end{bmatrix} = {{1} \over {\delta^2}} \begin{bmatrix} -2 & 1 & 0 & \cdots & 0 & 0 \\ 1 & -2 & 1 & \cdots & 0 & 0 \\ 0 & 1 & -2 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & -2 & 1 \\ 0 & 0 & 0 & \cdots & 1 & -2 \end{bmatrix} \begin{bmatrix} u (t, x_{1}) \\ u (t, x_{2}) \\ u (t, x_{3}) \\ \vdots \\ u (t, x_{m-2}) \\ u (t, x_{m-1}) \end{bmatrix} であり、ページを節約するために u′(t):=[ut(t,x1)ut(t,x2)ut(t,x3)⋮ut(t,xm−2)ut(t,xm−1)] \mathbf{u} ' (t) :=\begin{bmatrix} u_{t} (t, x_{1}) \\ u_{t} (t, x_{2}) \\ u_{t} (t, x_{3}) \\ \vdots \\ u_{t} (t, x_{m-2}) \\ u_{t} (t, x_{m-1}) \end{bmatrix}

Λ:=1δ2[−210⋯001−21⋯0001−2⋯00⋮⋮⋮⋱⋮⋮000⋯−21000⋯1−2] \Lambda := {{1} \over {\delta^2}} \begin{bmatrix} -2 & 1 & 0 & \cdots & 0 & 0 \\ 1 & -2 & 1 & \cdots & 0 & 0 \\ 0 & 1 & -2 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & -2 & 1 \\ 0 & 0 & 0 & \cdots & 1 & -2 \end{bmatrix}

u(t):=[u(t,x1)u(t,x2)u(t,x3)⋮u(t,xm−2)u(t,xm−1)] \mathbf{u} (t) := \begin{bmatrix} u (t, x_{1}) \\ u (t, x_{2}) \\ u (t, x_{3}) \\ \vdots \\ u (t, x_{m-2}) \\ u (t, x_{m-1}) \end{bmatrix} とし、u’(t)=Λu(t)\mathbf{u} ’ (t) = \Lambda \mathbf{u} (t)となる。今、時間ttのステップサイズをhhとし、tk=hkt_{k} = hkとして、いくつかの方法を適用して数値解を求めればよい。

■

教訓: スティッフな微分方程式

しかし、これはあくまで理論上の話で、実際に数値解を求めるには、かなりの困難が伴う。

Λ\Lambdaの固有値を知るために、Λ−λIm−1\Lambda - \lambda I_{m-1}の行列式を求めてみよう。

三重対角行列の行列式: 三重対角行列Xn:=[x10⋯001x1⋯0001x⋯00⋮⋮⋮⋱⋮⋮000⋯x1000⋯1x]X_{n} := \begin{bmatrix} x & 1 & 0 & \cdots & 0 & 0 \\ 1 & x & 1 & \cdots & 0 & 0 \\ 0 & 1 & x & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & x & 1 \\ 0 & 0 & 0 & \cdots & 1 & x \end{bmatrix}について、∣xn∣=Un(x2)\displaystyle | x_{n}| = U_{n} \left( {{x} \over {2}} \right)

Λ−λIm−1=1δ2[−2−δ2λ10⋯001−2−δ21⋯0001−2−δ2⋯00⋮⋮⋮⋱⋮⋮000⋯−2−δ21000⋯1−2−δ2] \Lambda - \lambda I_{m-1} = {{1} \over { \delta^2 }} \begin{bmatrix} -2 - \delta^2 \lambda & 1 & 0 & \cdots & 0 & 0 \\ 1 & -2 - \delta^2 & 1 & \cdots & 0 & 0 \\ 0 & 1 & -2 - \delta^2 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & -2 - \delta^2 & 1 \\ 0 & 0 & 0 & \cdots & 1 & -2 - \delta^2 \end{bmatrix} 三重対角行列である(m−1)×(m−1)(m - 1) \times (m-1)なので、行列式は det⁡(Λ−λIm−1)=1δ2Um−1(−2−δ2λ2) \det \left( \Lambda - \lambda I_{m-1} \right) = {{1} \over { \delta^2 }} U_{m-1} \left( {{-2 - \delta^2 \lambda } \over {2}} \right) のように求められる。ここで、Um−1U_{m-1}はm−1m-1次の第二種チェビシェフ多項式を意味する。

第二種チェビシェフ多項式の根と対称性:

  • [1]: Un(X)\displaystyle U_{n} (X)の根はxk=cos⁡(kn+1π)\displaystyle x_{k} = \cos \left( {{k} \over {n+1}} \pi \right)、k=1,⋯ ,nk=1, \cdots , n
  • [2]: Un(−x)=(−1)nUn(x)U_{n} (-x) = (-1)^{n} U_{n} (x)

括弧から−1-1を取り出せば (−1)m−1δ2Um−1(2+δ2λ2) {{(-1)^{m-1}} \over { \delta^2 }} U_{m-1} \left( {{2 + \delta^2 \lambda } \over {2}} \right) であり、固有値を求めるためには Um−1(2+δ2λ2) U_{m-1} \left( {{2 + \delta^2 \lambda } \over {2}} \right) だけ気にすれば十分である。Um−1(2+δ2λ2)\displaystyle U_{m-1} \left( {{2 + \delta^2 \lambda } \over {2}} \right)のjj番目の根はj=1,⋯ ,(m−1)j=1, \cdots , (m-1)に対して 2+δ2λj2=cos⁡(jmπ) {{2 + \delta^2 \lambda_{j}} \over {2}} = \cos \left( {{j} \over {m}} \pi \right) であり、両辺に22を掛けると 2+δ2λj=2cos⁡(jmπ)=2cos⁡(2⋅j2mπ)=2[1−2sin⁡2(jπ2m)] \begin{align*} \displaystyle 2+\delta^2 \lambda_{j} =& 2 \cos \left( {{j} \over {m}} \pi \right) \\ =& 2 \cos \left( 2 \cdot {{j} \over {2m}} \pi \right) \\ =& 2 \left[ 1 - 2 \sin^2 \left( {{ j \pi } \over { 2m }} \right) \right] \end{align*} したがって、固有値は λj=−4δ2sin⁡2(jπ2m) \lambda_{j} = - {{4} \over { \delta^2 }} \sin^2 \left( {{j \pi } \over {2m}} \right) のように求められる。しかし、 ∣λm−1∣=∣4δ2sin⁡2((m−1)π2m)∣≈4δ2=4m2 \begin{align*} \left| \lambda_{m-1} \right| =& \left| {{4} \over { \delta^2 }} \sin^2 \left( {{(m-1) \pi } \over {2m}} \right) \right| \\ \approx & {{4} \over {\delta^2}} \\ =& 4 m^2 \end{align*}

であるため、m=1δ\displaystyle m = {{1} \over {\delta}}が大きくなるほど、固有値の大きさも莫大になる。数値解をできるだけ正確に知るためにδ\deltaを小さくして[0,1][0,1]をたくさん分割した場合、計算は非常に難しくなる。たとえばm=100m=100にしてもu’(t)=Λu(t)\mathbf{u} ’ (t) = \Lambda \mathbf{u} (t)はスティッフstiffな問題となり、数値的に解くのはかなり困難である。

参照