2024-3S 計算数理演習（東京大学理学部・教養学部） [齊藤宣一] [top] [0] [1] [2] [3] [4] [5] [6] [7] [8] [UTOL]

8. 偏微分方程式：熱伝導，波動の数値計算

熱方程式

一様な材質でできた，長さが $L$ の針金の熱伝導現象は，次の偏微分方程式で記述される： \begin{equation} \tag{8.1a} \frac{\partial u}{\partial t}=\kappa \frac{\partial^2 u}{\partial x^2}\qquad (0 < x < L,~ t > 0). \end{equation} ここで，$u=u(x,t)$ は時刻 $t$ における点 $x$ での針金の温度である．また，$\kappa$ は針金の比熱や密度などから定まる正定数である．(8.1a)は，熱伝導方程式，熱拡散方程式，あるいは，単に，熱方程式と呼ばれる．

針金の温度変化を論じるには針金の端点（境界）での物理的な状況設定が必要である． もし，この両端が一定の温度，たとえば，$0$ に保たれているのであれば，$u$ は境界条件 \begin{equation} \tag{8.1b} u(0,t)=0,\quad u(L,t)=0\qquad (t > 0) \end{equation} を満たす．これを，Dirichlet境界条件と呼ぶ．

もし，端点で熱の出入りがない状態が保たれている，すなわち，断熱の状態であれば， 境界条件として， Neumann境界条件 \begin{equation*} \frac{\partial u}{\partial x}(0,t)=0,\quad \frac{\partial u}{\partial x}(L,t)=0 \qquad (t>0) \end{equation*} を採用しなけばならない．

温度変化の考察を $t=0$ から始めるとすると，この初期時刻での温度分布は既知であるとする必要がある．すなわち，$u$ には初期条件 \begin{equation} \tag{8.1c} u(x,0)=a\qquad (0\le x\le L) \end{equation} が課せられる．ただし，$a=a(x)$ は区間 $[0,L]$ 上で与えられた関数である．特に何も断らなければ，考えている境界条件を満たすような連続関数であるとしておく（実際には，この条件は，それほど本質的ではない）．

当面，熱方程式(8.1a)，Dirichlet境界条件(8.1b)，初期条件(8.1c)からなる，熱方程式の初期値境界値問題(8.1a,b,c)を扱う．Neumann境界条件については，次節で，詳しく扱う．

熱方程式の初期値境界値問題(8.1a,b,c)は，Fourierの方法を用いて，解を得ることができる．しかし，ここでは，数値的な方法による近似解法を考えたい．

差分近似

高等学校で学んだように，滑らかな関数 $v(x)$ の $x=\alpha$ における微分係数は， \[ v'(\alpha)=\frac{d v}{dx}(\alpha)=\lim_{h\to 0}\frac{v(\alpha+h)-v(\alpha)}{h} \] で定義される．したがって，$h$ が $0$ に十分近ければ， 分数 \begin{equation} \tag{8.2a} \frac{v(\alpha+h)-v(\alpha)}{h} \end{equation} の値は $v'(\alpha)$ に十分に近いことが期待される．実際， $h>0$ とし，$I$ を $\alpha$ と $\alpha\pm h$ を含むような有界な閉区間とする． このとき，Taylorの定理により， \[ v(\alpha+h)=v(\alpha)+v'(\alpha)h+\frac{1}{2}v''(\alpha+\theta h)h^2 \] を満たす $\theta\in [0,1]$ が存在するので，これを変形して， \begin{equation} \tag{8.2b} \left|v'(\alpha)-\frac{v(\alpha+h)-v(\alpha)}{h}\right|=\left|\frac{1}{2}v''(\alpha+\theta h)h\right| \le \frac12 h\max_{x \in I}|v''(x)| \end{equation} が得られる．これは，$v(x)$ が $I$ 上で $C^2$ 級であれば，$v'(\alpha)$と$[v(\alpha+h)-v(\alpha)]/h$ の差は $h$ に比例して小さくなることを示している．分数(8.2a)のことを，前進差分商，あるいは，前進Euler近似と呼ぶ．同じように，$v'(\alpha)$ を近似する量として， \begin{gather} \frac{v(\alpha)-v(\alpha-h)}{h},\tag{8.3}\\ \frac{v(\alpha+\frac{h}{2})-v(\alpha-\frac{h}{2})}{h},\quad \frac{v(\alpha+h)-v(\alpha-h)}{2h} \tag{8.4} \end{gather} を，それぞれ，後退差分商（後退Euler近似），中心差分商と呼ぶ．実際， \begin{align} \left|v'(\alpha)-\frac{v(\alpha)-v(\alpha-h)}{h}\right|%=\left|\frac{1}{2}v''(a-\theta h)h\right| & \le \frac12 h \max_{x \in I} |v''(x)|,\\ \left|v'(\alpha)-\frac{v(\alpha+\frac{h}{2})-v(\alpha-\frac{h}{2})}{h}\right| & \le \frac{1}{24}h^2 \max_{x \in I} |v^{(3)}(x)| \end{align} を導くことができる．

2階の微分係数 $v''(\alpha)$ の近似ついては， \begin{equation} \tag{8.5a} \frac{v(\alpha-h)-2v(\alpha)+v(\alpha+h)}{h^2} \end{equation} が得られる．これを，$v''(\alpha)$ の2階の中心差分商と呼ぶ．実際， \begin{equation} \tag{8.5b} \left|v''(\alpha)-\frac{v(\alpha-h)-2v(\alpha)+v(\alpha+h)}{h^2}\right| \le \frac{1}{12} h^2\max_{x\in I} |v^{(4)}(x)| \end{equation} を導くことができる．

熱方程式に対する陽的スキーム

あらためて，Dirichlet境界条件下での熱方程式の初期値境界値問題(8.1a,b,c)を考えよう．$N$ を正の整数として，$\displaystyle{h=\frac{L}{N+1}}$とおき，$x_i=ih$ $(i=0,1,\ldots,N+1)$ と定める．このとき，$x_0=0$，$x_{N+1}=L$であることに注意せよ．次に，正の定数$\tau$ を固定して，$t_n=n\tau$ $(n=0,1,\ldots)$とおく．この格子点の各点 $(x_i,t_n)$ で，初期値境界値問題の解 $u(x_i,t_n)$ の近似値 $u_i^n$ を求めたい．

$\frac{\partial u}{\partial t}(x_i,t_n)$ に前進差分商(8.2a)，$\frac{\partial^2 u}{\partial x^2}(x_i,t_n)$ に2階中心差分商(8.5a)を適用することで， (8.1a)に対する近似方程式として， \begin{equation} \frac{u_i^{n+1}-u_i^n}{\tau}=\kappa \frac{u_{i-1}^{n}-2u_i^{n}+u_{i+1}^n}{h^2}\qquad (1\le i\le N,\ n\ge 0) \tag{8.6a} \end{equation} が得られる．境界条件と初期条件は，単純に， \begin{gather} u_0^n=u_{N+1}^n=0 \qquad (n\ge 1), \tag{8.6b}\\ u_i^0=a(x_i) \qquad (0\le i\le N+1)\tag{8.6c} \end{gather} とする．

(8.6a,b)は， \begin{equation} \tag{8.7} \lambda\stackrel{\textrm{def.}}{=} \kappa \frac{\tau}{h^2} \end{equation} と置くことで， \begin{equation} \tag{8.8} \left\{ \begin{array}{l} u_0^{n+1}=0,\\ u_i^{n+1}=(1-2\lambda)u_i^n+\lambda (u_{i-1}^n+u_{i+1}^n) \qquad (1\le i\le N,n\ge 0),\\ u_{N+1}^{n+1}=0 \end{array} \right. \end{equation} と表現できる．したがって，近似方程式(8.6a,b,c)は，図8.2のような模式図（スキーム，scheme）で表現できる． すなわち．$u_i^{n+1}$ を計算する際には，「その真下の値」$u_i^n$ と「左右の下の値」$u_{i-1}^n,u_{i+1}^n$ を加重平均しているのである．このことから，(8.6a,b,c)のような差分近似から導かれた近似方程式を差分近似スキーム，あるいは，差分スキームと呼ぶ．特に，(8.6a,b,c)を陽的差分スキーム（explicit finite difference scheme），あるいは単に，陽的スキーム，陽的解法と呼ぶ． （この他に，「陰的」差分スキームがあるが，本講義では扱わない．）

注意. $h=L/(N+1)$ としているのは， $1\le i\le N$ に対して，$x_i$ を 区間 $[0,L]$ の内点にするためである． 結果的に，各 $t_n$ $(n\ge 1)$ に対して，$u_i^n$ $(1\le i\le N)$ が求めるべき未知数となる． $h=L/N$ とした場合には，$0=x_0 < x_1 < \cdots < x_{N-1} < x_{N}=L$ となり， $u_i^n$ $(1\le i\le N-1)$ が求めるべき未知数となる．実際のところ，陽的スキームを考える限りは，$h=L/(N+1)$ と定める恩恵は，さほどではない．あくまで，個人的な好みの問題であるので，$h=L/N$ としたい人は，（その定義に合わせて，しかるべく修正して）そのようにしても良い．

実際に，(8.6a.b.c)を計算してみよう．この節で使う予定のモジュールをすべてインポートしておく．

In [1]:  
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation, rc
from IPython.display import HTML
import plotly.graph_objects as go
import plotly.offline as pyo

(8.6a,b,c)を計算する関数を作成する． ただし，時間区間については，$0\le t\le T$ で考えるものとする． 次は，関数を実行する際に，ユーザが指定するものである：

• 方程式のパラメータ：$L=\mathrm{L}$，$T=\mathrm{T}$，$\kappa=\mathrm{kappa}$,
• 差分スキームのパラメータ：$N=\mathrm{N}$，$\lambda=\mathrm{lam}$,
• 可視化のため補助的に導入するパラメータ：$\mathrm{num}$.

$\mathrm{xgrid}=(x_0,x_1,\ldots,x_{N+1})\in\mathbb{R}^{N+2}$ とする． 「現在時刻」における $u_n^i$ を記憶する配列として， $\mathrm{u}=(u_0^n,u_1^n,\ldots,u_{N+1}^n)\in\mathbb{R}^{N+2}$ を導入する．すると， \[ \mathrm{u[1\mbox{:}N+1]}=(u_1^n,\ldots,u_{N}^n),\quad \mathrm{u[0\mbox{:}N]}=(u_0^n,\ldots,u_{N-1}^n),\quad \mathrm{u[2\mbox{:}N+2]}=(u_2^n,\ldots,u_{N+1}^n) \] であるから，(8.6a,b,c)，すなわち，(8.8)の更新式は， \[ \mathrm{u[1\mbox{:}N+1]}=(1-2\cdot \mathrm{lam})\cdot \mathrm{u[1\mbox{:}N+1]} +\mathrm{lam}\cdot(\mathrm{u[0\mbox{:}N]}+\mathrm{u[2\mbox{:}N+2]}) \] とすれば良い．ここで，右辺の $\mathrm{u[1\mbox{:}N+1]}$ などは， $(u_1^n,\ldots,u_{N}^n)$ などを， 左辺の $\mathrm{u[1\mbox{:}N+1]}$ は，$(u_1^{n+1},\ldots,u_{N}^{n+1})$ を意味することに注意せよ．

その上で，$\mathrm{tgrid}=(t_0,t_{\mathrm{step}},t_{2\mathrm{step}},\ldots,t_{\mathrm{num}\cdot\mathrm{step}})\in\mathbb{R}^{\mathrm{num}+1}$，および， \[ \mathrm{sol}= \begin{pmatrix} u_0^0 & u_1^0 & \cdots & u_N^0 & u_{N+1}^0 \\ u_0^{\mathrm{step}} & u_1^{\mathrm{step}} & \cdots & u_N^{\mathrm{step}} & u_{N+1}^{\mathrm{step}} \\ \vdots & \vdots & & \vdots & \vdots\\ u_0^{\mathrm{num}\cdot\mathrm{step}} & u_1^{\mathrm{num}\cdot\mathrm{step}} & \cdots & u_N^{\mathrm{num}\cdot\mathrm{step}} & u_{N+1}^{\mathrm{num}\cdot\mathrm{step}} \end{pmatrix} \in\mathbb{R}^{(\mathrm{num}+1)\times (N+2)} \] と定める．関数の計算結果としては， $\mathrm{xgrid}$，$\mathrm{tgrid}$，$\mathrm{sol}$ を返す．

以上を行う関数が，In [2]に示した heat_0dbc_ex1 である．

In [2]:  
#（斉次）熱方程式、零Dirichlet境界条件、陽的スキーム
def heat_0dbc_ex1(L, kappa, initial, N, T, lam, num):
    
  h = L/(N+1)
  tau = lam*h*h/kappa
  nmax = int(T/tau)
  step=int(max(1, nmax/num))
  
  xgrid=np.linspace(0.0, L, N+2)
  tgrid=np.array([0.0])
  u=np.vectorize(initial)(xgrid)
  sol=np.copy(u)
  
  for n in range(nmax):
    u[1:N+1] = (1-2*lam)*u[1:N+1]+lam*(u[0:N]+u[2:N+2])
    tnow = (n+1)*tau
    if (n+1)%step==0:
      sol=np.vstack((sol,u))
      tgrid=np.append(tgrid, tnow)

  return xgrid, tgrid, sol

次に可視化のために，関数を3つ用意する．

• In [3] plot_solution：図のプロット、一枚の図に、グラフを重ねる
• In [4] plot_solution3d： xtu空間での曲線の描画 （plotly.graph_objects モジュールを使う）
• In [5] plot_animation： アニメーションの作成（IPython.displayモジュールを使う．表示の際には，matplotlib.animationなどを使う）結果を自動的にgifに保存する際には，事前にGoogleドライブをマウントしておくこと．

In [3]:  
#図のプロット、一枚の図に、グラフを重ねる
def plot_solution(xgrid, uvect):
  # 描画に関するパラメータの設定
  umax = np.max(sol)
  umin = np.min(sol)
  ulength = umax - umin
  umax += 0.05 * ulength
  umin -= 0.05 * ulength
  # t = 0                  
  plt.plot(xgrid,uvect[0,:],color='#ff8c00')
  # t > 0
  for n in range(1,uvect.shape[0]):
    plt.plot(xgrid,uvect[n,:],color='#00bfff')
  
  plt.xlabel('x')
  plt.ylabel('u')
  plt.ylim(umin, umax)
  plt.grid('on')
  plt.show()

In [4]:  
# xtu空間での曲線の描画
# plotly.graph_objects，plotly.offline モジュールを使う
# xgrid.shape <= tgrid.shape でないと上手く描けない
def plot_solution3d(xgrid, tgrid, uvect):  
  # 3Dプロットの作成
  fig = go.Figure()
  
  Nt = tgrid.shape[0]
  # t = 0
  n = 0
  tval = tgrid[n]*np.ones(Nt)
  fig.add_trace(go.Scatter3d(x=xgrid, y=tval, z=uvect[n,:], mode='lines',line=dict(color='#ff8c00', width=2)))
  # t > 0 
  for n in range(1, Nt):
    tval = tgrid[n]*np.ones(Nt)
    fig.add_trace(go.Scatter3d(x=xgrid, y=tval, z=uvect[n,:], mode='lines',line=dict(color='#00bfff', width=2)))

  # レイアウトの設定
  fig.update_layout(scene=dict(xaxis_title='x', yaxis_title='time', zaxis_title='u'))
  fig.update_layout(scene=dict(camera=dict(eye=dict(x=1.5, y=-1.5, z=1))))
  fig.update_layout(showlegend=False)

  # プロットの表示
  pyo.plot(fig)
  with open('temp-plot.html', 'r') as file:
    html_content = file.read()
  display(HTML(html_content))
  # Google Colaboratoryでなければ、上の4行は、以下の1行で置き換えられる（はず）
  # fig.show()

In [5]: 
#アニメーションの作成
# IPython.displayモジュールを使う
def plot_animation(xgrid, tgrid, uvect):
  # 描画に関するパラメータの設定
  umax = np.max(sol)
  umin = np.min(sol)
  xmax = np.max(xgrid)
  xmin = np.min(xgrid)
  ulength = umax - umin
  utop = umax + 0.1*ulength
  umax += 0.05*ulength
  umin -= 0.05*ulength
  xmid = (xmax+xmin)/2
  xlength = xmax - xmin
  xmax += 0.05*xlength
  xmin -= 0.05*xlength
  
  # fig, axオブジェクトを作成
  fig, ax = plt.subplots()
  ims = []
  
  # t = 0
  n = 0
  im, = ax.plot(xgrid, uvect[n, :], color='#ff8c00')
  title = ax.text(xmid, utop, f'time={tgrid[n]:.4f}', ha='center', va='center', fontsize=15)
  ims.append([im, title])
    
  # t > 0
  for n in range(1, uvect.shape[0]):
    im, = ax.plot(xgrid, uvect[n, :], color='#00bfff')
    title = ax.text(xmid, utop, f'time={tgrid[n]:.4f}', ha='center', va='center', fontsize=15)
    ims.append([im, title])
    
  # 各軸のラベル
  ax.set_xlabel(r"$x$", fontsize=15)
  ax.set_ylabel(r"$u$", fontsize=15)
  # グラフの範囲を設定
  ax.set_xlim([xmin, xmax])
  ax.set_ylim([umin, umax])
  ax.grid(True)

  # ArtistAnimationにfigオブジェクトとimsを代入してアニメーションを作成
  return animation.ArtistAnimation(fig, ims)

これらを，次のIn [6]ように使う．熱方程式に対する差分法の最初の例としては，初期関数として， \[ a(x)=\min \{x ,\ L-x\} \] を用いるのが，慣例になっている． ここでもその慣例に従おう（この関数に対して，配列 $\mathrm{xgrid}$ を渡したときに，同じ次元の配列 $\mathrm{u}$ が返ってくるように，$\mathrm{u}=\mathrm{np.vectorize(initial)(xgrid)}$としている） ．

In [6]:  
# 初期値
def initial1(x):
  return min(x,1-x)
  
#パラメータの設定
L = 1.0
kappa = 1.0
N = 51
T = 1.0
lam = 0.4
num = 60

#熱方程式の計算
x, tn, sol = heat_0dbc_ex1(L, kappa, initial1, N, T, lam, num)

#図の表示
plot_solution(x, sol)

#3D図の表示
plot_solution3d(x, tn, sol)

#アニメーションを表示
anim = plot_animation(x, tn, sol)
#結果をheat.gifに保存（Googleドライブをマウントしておくこと！）
anim.save('/content/drive/MyDrive/Colab Notebooks/fig/heat.gif', writer='pillow')
rc('animation', html='jshtml')
plt.close()
anim

以下では，可視化にはplot_solutionのみを使う．

まずは，$L$，$\kappa$ を固定して，$N$ と $\lambda$ を色々変えて，計算をしてみよ．すると，例えば，$N=23$，$\lambda=0.51$，$T=0.18$として計算してみると，図8.4のようになる．

実は，(8.1)の解 $u$ は，初期関数 $a$ が連続で，$a(0)=a(L)=0$ 満たす限りにおいて， \begin{equation} \tag{8.9} \max_{0\le x\le L}|u(x,t)|\le \max_{0\le x\le L}|a(x)| \qquad (t\ge 0) \end{equation} を満たす（例えば，[藤田1977]の第1章を見よ）．これを，熱方程式の解の $L^\infty$ 安定性という． しかしながら，図8.4では，このような安定性は観察できない．すなわち，図8.4における振動現象は，近似の対象としている微分方程式の解の挙動を反映しているわけではない．

この例から，(8.9)と同様の性質は，陽的スキームの解については，一般には成り立たないことがわかる．それでは，どうすれば良いであろうか？ 実は， 対応する性質の成立を保証するためには，$\lambda$，すなわち，$h$，$\tau$ と $\kappa$ の間に， \begin{equation} \tag{8.10} \lambda\left(=\kappa\frac{\tau}{h^2}\right)\le \frac12 \end{equation} を仮定する必要がある．実際，(8.6a)の両辺の絶対値をとり，$1-2\lambda\ge 0$ に注意して， 三角不等式を使うと \begin{equation*} |u_{i}^{n+1}|\le (1-2\lambda)|u_{i}^n|+\lambda(|u_{i-1}^n|+|u_{i+1}^n|)\qquad (1\le i\le N) \end{equation*} となる．ここで，$M_n=\max_{1\le i\le N}|u_{i}^n|$ とおくと，この不等式に より，$|u_{i}^{n+1}|\le (1-2\lambda)M_n+\lambda(M_n+M_n)=M_n$ となるから， 両辺の $i$ についての最大値をとると，$M_{n+1}\le M_n$ を得る．したがって， $M_n\le M_0$，すなわち， \begin{equation} \tag{8.11} \max_{1\le i\le N}|u_{i}^n|\le \max_{1\le i\le N}|a(x_i)|\qquad (n\ge 0) \end{equation} が得られる．これを，差分解の $L^\infty$ 安定性， あるいは，離散 $L^\infty$ 安定性などという．特に，条件(8.10)が満たされている時には， 陽的スキームの解は発散しないことがわかる．

In [2]のheat_0dbc_ex1において， $N$ と $\tau$ の値でなく， $N$ と $\lambda$ の値をはじめに指定して，$\tau$ を，$\tau=\lambda h^2/\kappa$ で求めた理由は，(8.10)を容易に成立させるためである．

非斉次の熱方程式

針金に熱の湧き出しや吸収がある場合には，(8.1a)の代わりに，非斉次の熱方程式， \begin{equation} \tag{8.12} \frac{\partial u}{\partial t}=\kappa \frac{\partial^2 u}{\partial x^2}+f\qquad (0 < x < L,~ t > 0) \end{equation} を考えることになる．ここで，$f=f(x,t)$ は，与えられた関数である． 境界条件(8.1b)と初期条件(8.1c)は同じとする． このときには，(8.6a)を， \begin{equation} \tag{8.13} \frac{u_i^{n+1}-u_i^n}{\tau}=\kappa \frac{u_{i-1}^{n}-2u_i^{n}+u_{i+1}^n}{h^2}+ f(x_i,t_{n}) \qquad (1\le i\le N,\ n\ge 0), \end{equation} すなわち，(8.7)と同じ $\lambda$ を使って， \begin{equation*} u_i^{n+1} = (1-2\lambda)u_i^n+\lambda (u_{i-1}^{n}+u_{i+1}^n)+ \tau f(x_i,t_{n}) \qquad (1\le i\le N,\ n\ge 0) \end{equation*} と変えれば良い．(8.6b,c)はそのまま使える．

In [7]:  
#非斉次熱方程式、零Dirichlet境界条件、陽的スキーム
def heat_0dbc_ex2(L, kappa, initial, righthand, N, T, lam, num):
    
  h = L/(N+1)
  tau = lam*h*h/kappa
  nmax = int(T/tau)
  step=int(max(1, nmax/num))
  
  xgrid=np.linspace(0.0, L, N+2)
  tgrid=np.array([0.0])
  u=np.vectorize(initial)(xgrid)
  sol=np.copy(u)
  
  n=0
  tnow=n*tau
  for n in range(nmax):
    u[1:N+1] = (1-2*lam)*u[1:N+1]+lam*(u[0:N]+u[2:N+2]) + tau*righthand(xgrid[1:N+1], tnow)
    tnow = (n+1)*tau
    if (n+1)%step==0:
      sol=np.vstack((sol,u))
      tgrid=np.append(tgrid, tnow)

  return xgrid, tgrid, sol

$f$ と $a$ を， \begin{equation} \tag{8.14} f(x,t)=e^t(-x^4+x^3+12x^2-6x),\qquad a(x)=x^3(1-x) \end{equation} として計算をしてみよう．結果は，図8.5の通りである． 今度は，解 $u$ は，$0$ に減衰しない．

In [8]:  
# 初期値 a
def initial2(x):
  return x**3*(1.0 - x)

# 右辺の関数 f
def righthand2(x, t):
  return np.exp(t)*(-x**4 + x**3 + 12*x**2 - 6*x)
  
#パラメータの設定
L = 1.0
kappa = 1.0
N = 51
T = 1.0
lam = 0.4
num = 30

#熱方程式の計算
x, tn, sol = heat_0dbc_ex2(L, kappa, initial2, righthand2, N, T, lam, num)

#図の表示
plot_solution(x, sol)

次に，$f$ と $a$ を， \begin{equation} \tag{8.15} f(x,t)=f_3(x)=4\pi^2\sin(2\pi x),\qquad a(x)=x(1-x)\sin(4\pi x) \end{equation} として計算をしてみよう．

In [9]:  
# 初期値 a
def initial3(x):
  return np.sin(4.0*np.pi*x)*x*(1.0 - x)

# 右辺の関数 f
def righthand3(x, t):
  return 4.0*np.pi**2*np.sin(2.0*np.pi*x)
  
#パラメータの設定
L = 1.0
kappa = 1.0
N = 51
T = 1.0
lam = 0.4
num = 30

#熱方程式の計算
x, tn, sol = heat_0dbc_ex2(L, kappa, initial3, righthand3, N, T, lam, num)

#図の表示
plot_solution(x, sol)  

結果は，図8.6の通りである． このとき．数値解は，時間が十分に経過した際に，ある関数 $w=w(x)$ に漸近しているように見える．実際，（各点ごとに）$u(x,t)\to w(x)$ $(t\to\infty)$ を仮定すると，(8.12)と(8.1b)により，$w$ は， \[ 0=\kappa w''+f_3(x),\qquad w(0)=w(1)=0 \] を満たすはずである．したがって，$w(x)=\sin(2\pi x)$ であることがわかる．実際， $u$ が，$t\to\infty$ の際に，$w$ に $[0,1]$ で一様収束することを示すのは難しくない（例えば，[藤田1977]を見よ）．

誤差の考察

非斉次熱方程式の初期値境界値問題(8.12)，(8.1b,c)と，その陽的スキーム(8.13)，(8.6b,c)の誤差について考察する．その前に，まず，(8.13)，(8.6b,c)の解 $u_i^n$ についても， (8.11)に対応した結果を導いておこう．実際， $\lambda\le 1/2$ である限りにおいて， \begin{equation*} |u_{i}^{n+1}|\le (1-2\lambda)|u_{i}^n|+\lambda(|u_{i-1}^n|+|u_{i+1}^n|)+ \tau |f(x_i,t_n)| \quad (1\le i\le N) \end{equation*} と評価できるので，離散 $L^\infty$ 安定性 \begin{equation} \tag{8.16} \max_{1\le i\le N}|u_{i}^n|\le \max_{1\le i\le N}|a(x_i)| +\tau \sum_{k=0}^{n-1} \max_{1\le i\le N} |f(x_i,t_k)| \qquad (n\ge 1) \end{equation} が得られる．これを見やすくするために，一般に，$\boldsymbol{v}=(v_1,\ldots,v_N)\in\mathbb{R}^{N}$ に対して，$\infty$ ノルム \[ \|\boldsymbol{v}\|_\infty = \max_{1\le i\le N}|v_{i}| \] を導入し，$\boldsymbol{u}^{(n)}=(u_1^n,\ldots,u_N^n)$， $\boldsymbol{a}=(a(x_1),\ldots,a(x_N))$， $\boldsymbol{f}^{(n)}=(f(x_1,t_n),\ldots,f(x_i,t_n))$ とおく．すると，(8.16)は， \begin{equation} \tag{8.16} \|\boldsymbol{u}^{(n)}\|_\infty\le \|\boldsymbol{a}\|_\infty +\tau \sum_{k=0}^{n-1} \|\boldsymbol{f}^{(k)}\|_\infty \qquad (n\ge 1) \end{equation} と書ける．

さて，改めて，誤差 $e_i^n = u_i^n - u(x_i,t_n)$ について考察しよう． 常微分方程式の誤差解析の時と同様に，この誤差が，どのような差分スキームを満たすのかを考えることが出発点となる．$U_i^n=u(x_i,t_n)$ とおくと， \begin{align*} &\frac{e_i^{n+1}-e_i^n}{\tau}-\kappa \frac{e_{i-1}^{n}-2e_i^{n}+e_{i+1}^n}{h^2}\\ &=\left(\frac{u_i^{n+1}-u_i^n}{\tau}-\kappa \frac{u_{i-1}^{n}-2u_i^{n}+u_{i+1}^n}{h^2}\right) -\left(\frac{U_i^{n+1}-U_i^n}{\tau}-\kappa \frac{U_{i-1}^{n}-2U_i^{n}+U_{i+1}^n}{h^2}\right)\\ &= \underbrace{f(x_i,t_n)- \left(\frac{U_i^{n+1}-U_i^n}{\tau}-\kappa \frac{U_{i-1}^{n}-2U_i^{n}+U_{i+1}^n}{h^2}\right)}_{=r_i^n\mbox{とおく}}\\ &=r_i^n. \end{align*} すなわち，$e_i^n$ は， \[ \frac{e_i^{n+1}-e_i^n}{\tau}=\kappa \frac{e_{i-1}^{n}-2e_i^{n}+e_{i+1}^n}{h^2}+ r_i^n \qquad (1\le i\le N,\ n\ge 0) \] と境界条件 $e_0^n=e_{N+1}^n=0$，および，初期条件 $e_i^0 = 0$ を満たす． したがって，$\lambda$ が(8.10)を満たすならば，離散 $L^\infty$ 安定性(8.16)により， \[ \|\boldsymbol{e}^{(n)}\|_\infty\le \tau \sum_{k=0}^{n-1} \|\boldsymbol{r}^{(k)}\|_\infty \qquad (n\ge 1) \] が得られる．ただし，$\boldsymbol{e}^{(n)}=(e_1^n,\ldots,e_N^n)$，$\boldsymbol{r}^{(n)}=(r_1^n,\ldots,r_N^n)$ である．

一方で，$r_i^n$ は，方程式(8.12)を使って， \begin{align*} r_i^n &= \frac{\partial u}{\partial t}(x_i,t_n)-\kappa \frac{\partial^2 u}{\partial x^2}(x_i,t_n)-\left(\frac{U_i^{n+1}-U_i^n}{\tau}-\kappa \frac{U_{i-1}^{n}-2U_i^{n}+U_{i+1}^n}{h^2}\right)\\ &= \left( \frac{\partial u}{\partial t}(x_i,t_n)-\frac{U_i^{n+1}-U_i^n}{\tau} \right)-\kappa \left(\frac{\partial^2 u}{\partial x^2}(x_i,t_n) -\frac{U_{i-1}^{n}-2U_i^{n}+U_{i+1}^n}{h^2}\right) \end{align*} と書ける．ここで，$T>0$ を任意に固定して，$Q=[0,L]\times [0,T]$ とおく．すると，差分商に対する誤差評価(8.2b)と(8.5b)を用いて， \[ |r_i^n|\le \frac{1}{2}\tau \max_{(x,t)\in Q} \left|\frac{\partial^2 u}{\partial t^2}(x,t)\right| +\frac{\kappa}{12}h^2 \max_{(x,t)\in Q} \left|\frac{\partial^4 u}{\partial x^4}(x,t)\right| \] と評価できる．したがって，さらに， \[ \mathcal{M}= \frac{1}{2} \max_{(x,t)\in Q} \left|\frac{\partial^2 u}{\partial t^2}(x,t)\right| +\frac{\kappa}{12}\max_{(x,t)\in Q} \left|\frac{\partial^4 u}{\partial x^4}(x,t)\right| \] とおくと， \[ \|\boldsymbol{r}^{(n)}\|_\infty \le (\tau+h^2) \mathcal{M} \qquad (0\le t_n\le T) \] が成り立つ．これに，上で導いた， $\|\boldsymbol{e}^{(n)}\|_\infty$ と$\|\boldsymbol{r}^{(n)}\|_\infty$ の関係式を使うと，次の定理が証明できたことになる．

定理． $T>0$ を任意に固定する． 非斉次熱方程式の初期値境界値問題(8.12)，(8.1b,c)の解 $u$ は，$\mathcal{M}<\infty$ となるほど滑らかであるとする．$\lambda$ について，(8.10)を仮定する． このとき，$u$ と陽的スキーム(8.13)，(8.6b,c)の解 $u_i^n$ の誤差は， \[ \max_{0\le t_n\le T}\|\boldsymbol{e}^{(n)}\|_\infty\le T\mathcal{M}(\tau+h^2) \] を満たす．

上の計算を実験的に確かめてみよう． そのために，(8.14)の $f$ と $a$ を使う．このとき，厳密解は， \[ u(x,t)=e^tx^3(1-x) \] となる．$\tau=\lambda h^2/\kappa$ とするため，誤差の挙動としては，$O(h^2)$ が予想される．

In [10]:  
#誤差の計算、非斉次熱方程式、零Dirichlet境界条件、陽的スキーム
def error_heat(L, kappa, initial, righthand, exact, T, lam):
    
  N = 10
  num = 200

  kmax = 6
  hv = np.zeros(kmax)
  ev = np.zeros(kmax)

  for k in range(kmax): 
    N = N*2
    x, tn, sol = heat_0dbc_ex2(L, kappa, initial, righthand, N, T, lam, num)
    err = 0.0
    for n in range(tn.shape[0]):
      tval = tn[n]
      err = max(err, np.linalg.norm(sol[n,:] - exact(x, tval), ord=np.inf))
    hv[k] = L/(N+1)
    ev[k] = err

  rate=(np.log(ev[1:]) - np.log(ev[:-1])) / (np.log(hv[1:]) - np.log(hv[:-1]))
  
  return hv, ev, rate

In [11]:  
# 初期値 a
def initial2(x):
  return x**3*(1.0 - x)

# 右辺の関数 f
def righthand2(x, t):
  return np.exp(t)*(-x**4 + x**3 + 12*x**2 - 6*x)
  
# 厳密解 u
def exact2(x, t):
  return np.exp(t)*(1.0 - x)*x**3
  

#パラメータの設定
L = 1.0
kappa = 1.0
N = 51
T = 1.0
lam = 0.4
num = 30

#熱方程式の計算
hv, ev, rate = error_heat(L, kappa, initial2, righthand2, exact2, T, lam)

#収束の速さの出力
for i in range(rate.shape[0]-1):
  print(f'{hv[i+1]:.3f}, {rate[i]:.3f}')

#結果の描画（両対数目盛）
plt.plot(hv, ev, 'bo-')   
plt.xscale('log')
plt.yscale('log')
plt.legend(['explicit'])  
plt.xlabel('h')
plt.ylabel('error')
plt.grid('on')
plt.gca().set_aspect('equal', adjustable='box') 
plt.show()
  
Out [11]: 
0.024, 1.995
0.012, 1.999
0.006, 2.000
0.003, 2.000

波動方程式

ギターの弦のように真っ直ぐに張った弦の横振動は，波動方程式 \begin{align} & \frac{\partial^2 u}{\partial t^2}=c^2 \frac{\partial^{2} u}{\partial x^{2}} && (0 < x < L, t>0),\tag{8.17a}\\ &u(0,t)=u(L,t)=0 && ( t > 0),\tag{8.17b}\\ &u(x,0)=f(x),\quad \frac{\partial u}{\partial t}(x,0)=g(x) && (0\le x\le L)\tag{8.17c} \end{align} で記述される．$u=u(x,t)$ は弦の鉛直方向の変位，$c>0$ は伝搬速度， $f$ は弦の初期形状，$g$ は弦の初期速度である．

引き続き，格子点の集合 $(x_i,t_n)$ $(0\le i\le N+1,~ n\ge 0)$ を考える． ただし，$N$ を正の整数として，$h=L/(N+1)$ とおき，$x_i=ih$ $(i=0,1,\ldots,N+1)$ と，正の定数 $\tau$ を固定して，$t_n=n\tau$ $(n=0,1,\ldots)$ としている．

熱方程式の場合と同様に考えると，(8.17a)の差分近似は， \begin{equation} \tag{8.18a} \frac{u_{i}^{n-1}-2u_{i}^n+u_{i}^{n+1}}{\tau^2} =c^2\frac{u_{i-1}^n-2u_{i}^n+u_{i+1}^n}{h^2}\qquad (1\le i\le N,\ n\ge 1) \end{equation} とすればよい．境界条件と一つ目の初期条件は， \begin{gather} u_{0}^n=u_{N+1}^n=0\qquad (n\ge 1),\tag{8.18b}\\ u_{i}^0=f(x_i) \qquad (0\le i\le N+1)\tag{8.18c} \end{gather} とする．二つ目の初期条件は前進差分商を応用して， \begin{equation} \tag{8.19a} \frac{u_{i}^1-u_{i}^0}{\tau}=g(x_i) \qquad (0\le i\le N+1) \end{equation} とするのが簡単であるが，次のように考えてもよい．すなわち，$t=0$ を中心にして，中心差分商を適用して，（仮に $n=-1$ を認めて） \begin{equation} \tag{8.19b} \frac{u_{i}^1-u_{i}^{-1}}{2\tau}=g(x_i) \qquad (0\le i\le N+1) \end{equation} が得られる．ここに現れる $u_{i}^{-1}$ を消去するために，$n=0$ において差分方程式(8.18a) の成立を要請すると， \begin{equation*} \frac{u_{i}^1-2u_{i}^0+u_{i}^{-1}}{\tau^2} =c^2\frac{u_{i-1}^0-2u_{i}^0+u_{i+1}^0}{h^2} \end{equation*} となる．ここで， \begin{equation} \tag{8.20} \lambda \stackrel{\textrm{def.}}{=} c\frac{\tau}{h} \end{equation} と定義する．熱方程式の(8.7)とは異なるので注意せよ． そして，上の関係式から，$u_{i}^{-1}$ を消去すると， \begin{align*} u_{i}^1 &=2u_{i}^0-u_{i}^{-1}+\lambda^2(u_{i-1}^0-2u_{i}^0+u_{i+1}^0)\\ &=2u_{i}^0-u_{i}^1+2g(x_i)\tau +\lambda^2(u_{i-1}^0-2u_{i}^0+u_{i+1}^0)\\ &=2(1-\lambda^2)f(x_i)-u_{i}^1+2g(x_i)\tau +\lambda^2(f(x_{i-1})+f(x_{i+1})), \end{align*} すなわち，二つ目の初期条件の近似として， \begin{equation} \tag{8.18d} u_{i}^1 =g(x_i)\tau +(1-\lambda^2)f(x_i)+\frac{\lambda^2}{2}(f(x_{i-1})+f(x_{i+1})) \qquad (0\le i\le N+1) \end{equation} が得られる．

まとめると， 波動方程式に対する初期値境界値問題(8.17a,b,c)に対する差分解法は，次のようになる．まず， (8.18c)と(8.18d)で， $\{u_{i}^0\}_{1\le i\le N}$ と $\{u_{i}^1\}_{1\le i\le N}$ を計算し，次に，$n\ge 1$ に対して， \begin{equation*} u_{i}^{n+1} =2(1-\lambda^2)u_{i}^n+\lambda^2(u_{i-1}^n+u_{i+1}^n)-u_{i}^{n-1}\qquad (1\le i\le N) \end{equation*} で $\{u_{i}^{n+1}\}_{1\le i\le N}$ を計算するわけである．

In [12]:  
#非斉次波動程式、零Dirichlet境界条件、陽的スキーム
def wave_0dbc_ex(L, c, funcf, funcg, N, T, lam, num):
    
  h = L/(N+1)
  tau = lam*h/c
  nmax = int(T/tau)
  step=int(max(1, nmax/num))
  
  xgrid=np.linspace(0.0, L, N+2)
  tgrid=np.array([0.0])
  unew = np.zeros_like(xgrid)
  
  # n = 0
  u = np.vectorize(funcf)(xgrid)
  n = 0
  tnow = n*tau
  sol = np.copy(u)
  upast = np.copy(u)
  
  # n = 1
  u[1:N+1]=tau*np.vectorize(funcg)(xgrid[1:N+1])+(1-lam**2)*u[1:N+1]+(lam**2/2.0)*(u[0:N]+u[2:N+2])
  n = 1 
  tnow=n*tau
  sol=np.vstack((sol,u))
  tgrid=np.append(tgrid, tnow)
  
  # n >= 2 
  for n in range(1, nmax, 1):
    unew[1:N+1] = 2*(1-lam**2)*u[1:N+1] + lam**2*(u[0:N]+u[2:N+2]) - upast[1:N+1]
    upast=np.copy(u)
    u=np.copy(unew)  
    tnow=(n+1)*tau 
    if (n+1)%step==0:
      sol=np.vstack((sol,u))
      tgrid=np.append(tgrid, tnow)

  return xgrid, tgrid, sol

In [13]:  
# 初期値 f
def displacement(x):
  return np.sin(2.0*np.pi*x)
      
# 初期値 g
def velocity(x):
  return 0.0
  
#パラメータの設定
L = 1.0
c = 1.0
N = 51
T = 2
lam = 1.0
num = 60

#波動方程式の計算
x, tn, sol = wave_0dbc_ex(L, c, displacement, velocity, N, T, lam, num)

#図の表示
plot_solution(x, sol)

#アニメーションを表示
anim = plot_animation(x, tn, sol)
#結果をheat.gifに保存（Googleドライブをマウントしておくこと！）
anim.save('/content/drive/MyDrive/Colab Notebooks/fig/wave.gif', writer='pillow')
rc('animation', html='jshtml')
plt.close()
anim

$f$ と $g$ を， \begin{equation} \tag{8.21} f(x,t)=\sin (2\pi x),\qquad g(x,t)=0 \end{equation} とした際の計算結果を．図8.8a, bに示す． $c=1$，$N=51$，$\lambda=1$，$T=2$ としている．

しかしながら，図8.9に示した通り，$\lambda=1.1$ とすると， 差分スキームの解は破綻してしまう．

実は，$\lambda$ が， \begin{equation} \tag{8.22} \lambda\le 1 \end{equation} を満たすならば，差分スキームの解は（ある意味において)安定であり，また，波動方程式の解に収束する．この(8.22)を，Courant--Friedrichs--Lewy (CFL) 条件という（これは，[CFL]を起源としている）．

しかしながら，熱方程式の(8.16)のような顕著な安定性が得られないため，誤差解析は，案外厄介である．興味のある人は，[齊藤2023]や[John]を見よ．したがって，本講義では，数値実験のみで誤差の挙動を確かめておく．$c=1$ として，(8.21)の $f$ と $g$ を用いたとき， \[ u(x,t)=\cos (2\pi t)\sin(2\pi x) \] が厳密解となる． Out [15]により，誤差の挙動は $O(\tau^2+h^2)$ であることが予想される（そして実際にこの予想は正しい）．

In [14]:  
#誤差の計算、波動方程式、零Dirichlet境界条件、陽的スキーム
def error_wave(L, c, funcf, funcg, exact, T, lam):
    
  N = 10
  num = 100

  kmax = 6
  hv = np.zeros(kmax)
  ev = np.zeros(kmax)

  for k in range(kmax): 
    N = N*2
    x, tn, sol = wave_0dbc_ex(L, c, funcf, funcg, N, T, lam, num)
    err = 0.0
    for n in range(tn.shape[0]):
      tval = tn[n]
      err = max(err, np.linalg.norm(sol[n,:] - exact(x, tval), ord=np.inf))
    hv[k] = L/(N+1)
    ev[k] = err
 
  rate=(np.log(ev[1:]) - np.log(ev[:-1])) / (np.log(hv[1:]) - np.log(hv[:-1]))

  return hv, ev, rate

In [15]:  
# 初期値 f
def displacement(x):
  return np.sin(2.0*np.pi*x)
      
# 初期値 g
def velocity(x):
  return 0.0
  
# 厳密解 u
def exact_wave(x, t):
  return np.sin(2.0*np.pi*x)*np.cos(2.0*np.pi*t)
  
#パラメータの設定
L = 1.0
c = 1.0
T = 1.0
lam = 0.9

#波動方程式の計算
hv, ev, rate = error_wave(L, c, displacement, velocity, exact_wave, T, lam)

#収束の速さの出力
for i in range(rate.shape[0]-1):
  print(f'{hv[i+1]:.3f}, {rate[i]:.3f}')

#結果の描画（両対数目盛）
plt.plot(hv, ev, 'bo-')   
plt.xscale('log')
plt.yscale('log')
plt.legend(['explicit'])  
plt.xlabel('h')
plt.ylabel('error')
plt.grid('on')
plt.gca().set_aspect('equal', adjustable='box') 
plt.show()  

Out [15]: 
0.024, 2.004
0.012, 2.000
0.006, 2.000
0.003, 2.000

注意. In [15] において，$\lambda=1$ とすると，「精度が良すぎて」逆に精度が観察できない．各自，試してみよ．

注意. (8.22)を，Courant--Friedrichs--Lewy (CFL) 条件と言うと述べた． 一方で，熱方程式に対する(8.10)も，CFL条件と言われることがあるが，私個人は，適切な用語法とは思えないため，使わない．

問題

1. 上記の入力と出力をすべて確かめよ．
2. パラメータ値を変えて，計算を再試行せよ．
3. 例題を自分で設定して，上記の方法を試せ．

課題

課題8.1

熱方程式の初期値境界値問題(8.1a,b,c)に対して，厳密解を自分で作り，$\lambda=1/6$ のとき， 陽的スキーム(8.6a,b,c)の誤差について調べよ． 可能ならば，観察した事実を，数学的に説明せよ． さらに， 非斉次の場合(8.12)と(8.13)についても同様の観察をせよ．

課題8.2

熱方程式(8.1a)に対して，初期条件(8.1c)と，非斉次のDirichlet境界条件 \[ u(0,t)=1,\quad u(L,t)=0\qquad (t>0) \] を与えた時の初期値境界値問題を考え，差分スキームを導け．また， $t_n\to \infty$ の際の数値解 $u_i^n$ の挙動を観察せよ．

課題8.3

初期条件 $\frac{\partial u}{\partial t}(x,0)=g(x)$ の近似を，(8.19b)でなく(8.19a)とした場合の，波動方程式の初期値境界値問題(8.17a,b,c)に対する差分スキームを導け．また，実際に計算を遂行して，誤差の挙動を調べよ．

課題8.4

波動方程式の初期値境界値問題(8.17a,b,c)において，境界条件を，周期境界条件 \[ u(0,t)=u(L,t),\quad \frac{\partial u}{\partial x}(0,t)= \frac{\partial u}{\partial x}(L,t)\qquad (t>0) \] に代えた問題に対して差分法を適用せよ（$f$ と $g$ は境界条件に応じて，適切に選ぶこと）．

ヒント: この場合は，$x_i=ih$ $(i=0,1,\ldots,N)$，$h=L/N$ と定めた方が良い． そして，$u_i^n\approx u(x_i,t_n)$ に対して，境界条件により，$u_i^n=u_{i+N}^n$ $(\forall i)$ とする．

参考文献

• [齊藤2023] 齊藤宣一，偏微分方程式の計算数理，共立出版，2023年
• [藤田1977] 藤田宏, 池部晃生, 犬井鉄郎, 高見穎郎， 数理物理に現れる偏微分方程式I, 岩波講座基礎数学, 岩波書店, 1977年（岩波オンデマンドブックス版，2019年）
• [John] F. John, Lectures on advanced numerical analysis, Gordon and Breach Science Publishers, New York-London-Paris, 1967. 【和訳】藤田宏[訳]，名取亮[訳]，数値解析講義，産業図書，1975年．
• [CFL] R. Courant, K. Friedrichs, and H. Lewy, On the partial difference equations of mathematical physics, IBM J. Res. Develop., Vol. 11, pp. 215--234, 1967. 【原著】Über die partiellen Differenzengleichungen der mathematischen Physik. Math. Ann. 100 (1928), no. 1, 32--74.