2025-3A 数理情報学（東京大学教養学部） [齊藤宣一] [top] [A] [B] [C] [UTOL]

C. 非線形保存則とHamilton-Jacobi方程式に対する差分スキーム

非線形保存則に対する差分スキーム

非線形保存則の差分解法を説明する．講義では， \begin{align*} u_t+\frac{\partial}{\partial x}f(u)&=0, && (x,t)\in \mathbb{T}\times (0,T),\\ u(x,0)&={u}_0(x), && x\in \mathbb{T} \end{align*} を扱った．ただし，$\mathbb{T}=\mathbb{R}/\mathbb{Z}$，$\partial_x=\partial/\partial x$であり， $f:\mathbb{R}\to\mathbb{R}$， ${u}_0:\mathbb{R}\to\mathbb{R}$ は与えられた関数である．特に，$f$ は $C^1$ 級の関数とする．計算実習では，これを一般化して次の初期値問題を計算することにしよう： ：\begin{align} u_t+\frac{\partial}{\partial x}f(u)&=g(x,t), && (x,t)\in \mathbb{T}\times (0,T),\tag{C.1a}\\ u(x,0)&={u}_0(x), && x\in \mathbb{T}.\tag{C.1b} \end{align} ここで，$g(x,t)$ は与えられた連続関数である．

$0 < N,m\in\mathbb{Z}$ を取り，$h=1/N$，$\tau=T/M$， $x_i=ih$，$t_n=n\tau$ とする．$u^n_i\approx u(x_i,t_n)$ が求めるべき近似解である． 以下，何も断らなくても， \[ u_{i+N}^n=u_{i}^n \] とする．さらに， \begin{equation*} \lambda \stackrel{\textrm{def.}}{=}\frac{\tau}{h} \end{equation*} とおく．

(C.1)に対するLax–Friedrichs (LF) 差分スキームは， \begin{equation*} \frac{u_i^{n+1}-\frac{u^n_{i+1}+u^n_{i-1}}{2}}{\tau} +\frac{f(u^n_{i+1})-f(u^n_{i-1})}{2h}=g(x_i,t_n) \end{equation*} である．すなわち， \begin{align} u_i^{n+1} &= u^n_i-\lambda \left[ \frac{f(u^n_{i+1})-f(u^n_{i-1})}{2} -\frac{1}{2\lambda}(u^n_{i-1}-2u^n_{i}+u^n_{i+1}) \right]+\tau g(x_i,t_n) \nonumber \\ &= u_i^n-\lambda\Delta^+ \left[ \frac{f(u_i^n)+f(u_{i-1}^n)}{2} -\frac{1}{2\lambda}(u_i^n-u_{i-1}^n)\right]+\tau g(x_i,t_n) \tag{C.2} \end{align} である．なお，$\Delta^{+} v_i=v_{i+1}-v_i$ と表記している．

そして，初期条件は， \begin{equation*} u^{0}_i={u}_0(x_i) \end{equation*} とする（このことは，以後，いちいち断らない）．

風上差分スキームを述べるために，次を仮定する： \begin{equation} \exists \gamma\in\mathbb{R},\quad f'(\gamma)=0, \quad f'(s) \begin{cases} \ge 0 & (s > \gamma),\\ \le 0 & (s < \gamma). \end{cases} \end{equation} そして， \begin{equation} \label{eq:ncl-upwind-theta} \theta(s)= \begin{cases} 0 & (s>\gamma),\\ 1 & (s\le \gamma) \end{cases} \end{equation} とおき，混乱の恐れがない限り， \begin{equation*} \theta_i=\theta_i^n=\theta(u_i^n) \end{equation*} とかく．このとき，風上差分スキームは， \begin{equation} \label{eq:ncl-upwind2} u_i^{n+1}=u_i^n-\lambda \left[ \begin{array}{l} \theta_{i+1} (f(u_{i+1}^n)-f(\gamma))-\theta_i (f(u_i^n)-f(\gamma))+\\ (1-\theta_{i})(f(u_{i}^n)-f(\gamma))-(1-\theta_{i-1})(f(u_{i-1}^n)-f(\gamma)) \end{array} \right]+\tau g(x_i,t_n) \end{equation} となるのであった（詳細は講義で説明した）．

用語の確認をしておこう．$g\equiv 0$ とする．次の形で表現できる差分スキームを，保存型の差分スキームと言う： \begin{align} u_i^{n+1} &=u_i^n-\lambda\Delta^+ \varphi (u_{i-1}^n,u_i^n) \\ &=u_i^n-\lambda[\varphi (u_{i}^n,u_{i+1}^n)-\varphi (u_{i-1}^n,u_i^n)]. \end{align} この関数 $\varphi$ を数値流束と言う． LF差分スキームと風上差分スキームは保存型の差分スキームである．さらに，保存型の差分スキームを \[ u_i^{n+1}=\Phi(u_{i-1}^n,u_i^n,u_{i+1}^n) \] の形に書いたとき，それが区間 $[a,b]$ において単調であるとは， $u_{i-1},u_i,u_{i+1}\in [a,b]$ の範囲にある $u_{i-1},u_i,u_{i+1}$ について，関数 $\Phi(u_{i-1},u_i,u_{i+1})$ が各変数について非減少関数になることである（他の2つの変数を固定したときにある変数について非減少ということ）． 単調な保存型の差分スキームが，単調性，安定性，順序保存性などとても良い性質を持つことは講義で説明した． LF差分スキームと風上差分スキームは， Courant-Friedrichs-Lewy (CFL)条件 \begin{equation} \tag{C.3} \lambda \max_{s\in [a,b]}|f'(s)|\le 1 \end{equation} の下で，$[a,b]$ において単調になることを講義で証明した．

利用予定のモジュールを最初にまとめてインポートしておく．

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

In [2]:  
from google.colab import drive
drive.mount('/content/drive')

可視化のためのプログラムは，前回までと同じものを使う．

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]:  
#アニメーションの作成
# 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)

LF差分スキームを実際に計算してみよう．それには，線形移流方程式と同様に考えれば良い．

しかし，CFL条件(C.3)の成立には気を配る必要がある．まず，$a,b$ は，$a\le u_0(x)\le b$ を満たすようにとることにしよう．

これを実現するために，次の様に考える．まず，初期関数 $u_0(x)$ の最小値 $s_{\min}$ と最大値 $s_{\max}$ を計算する．その上で， $\Delta s= (s_{\max}-s_{\min})/1000$，$s_j= s_{\min}+j\Delta s$ $(j=0,\ldots,1000)$ とし， \[ \left|\frac{f(s_{j+1})-f(s_j)}{\Delta s}\right|\qquad (j=0,\ldots,1000-1) \] の最大値 $\mathrm{d}f_{\max}$ を取得する．これを計算するのがIn [5]であり，Out [6] のように計算できる．

In [5]:  
def find_maximum(ff, init):
  s0 = np.min(init)
  s1 = np.max(init)
  Ns = 1000
  xs = np.linspace(s0, s1, Ns)
  ys = ff(xs)
  dfs = (ys[1:] - ys[:-1])*Ns/(s1-s0)
  return s0, s1, np.max(np.abs(dfs))

In [6]:
def func_f0(s):
  return -(s-1)**2 + 3

def func_init(x):
  return np.cos(x)

x0 = 0
x1 = np.pi
xj = np.linspace(x0, x1, 100)
yj = func_init(xj)

s0, s1, ff_max = find_maximum(func_f0, yj)
print(s0, s1, ff_max)

Out [6]: 
-1.0 1.0 4.0019999979958865

(C.3)を計算するプログラムでは，$\mathrm{d}f_{\max}$ を $\max_{s\in [a,b]}|f'(s)|$ の代用として扱う． はじめに，$\mu=\lambda \max_{s\in [a,b]}|f'(s)|$ の値を，$0 < \mu\le 1$ の範囲で指定しておく． そして，$h$ を $h=1/N$ で定めた後に，$\tau=\mu h /\mathrm{d}f_{\max}$ および $\lambda=\tau/h$ とする．それが，In [7]である．

In [7]:  
def conservationlaw_cf(flux, initial, righthand, N, T, mu, num):
  # x-grids
  L = 1.0
  h = L/N
  xgrid = np.linspace(0.0, L, N)
  xtmp = xgrid[:-1].copy()

  # initial setteings
  tgrid = np.array([0.0])
  u = initial(xtmp)
  sol = np.copy(np.append(u,u[0]))
  
  # CFL condition
  smin, smax, abs_df = find_maximum(flux, u)

  # t-grids
  tau = mu*h/abs_df
  lam = tau/h
  nmax = int(T/tau)
  step = int(max(1, nmax/num))

  # iteration 
  const = 0.5/lam
  n = 0
  tnow = n*tau
  for n in range(nmax):
    uflux = 0.5*(flux(u) + flux(np.roll(u,1))) - const * (u - np.roll(u,1)) 
    u = u - lam*(np.roll(uflux,-1) - uflux) + tau*righthand(xtmp, tnow)
    tnow = (n+1)*tau
    if (n+1)%step==0:
      sol=np.vstack((sol,np.append(u,u[0])))
      tgrid=np.append(tgrid, tnow)

  return xgrid, tgrid, sol

例題を設定する．まず， \[ f(s)=\frac{1}{2}s^2 \] とする．その上で， 厳密な解が \[ u(x,t)=\sin(2\pi (x-t))+1.2 \] となるように，$u_0(x)$ と $g(x,t)$ を定める．In [8]のようにこれを実行する．

In [8]: 
# 初期値
def initial1(x):
  return np.sin(2*np.pi*x) + 1.2

# 流束関数
def flux1(s):
  return 0.5*s*s

#右辺の関数
def righthand1(x,t):
  g0 = 1.2
  g1 = np.sin(2*np.pi*(x-t)) + g0
  gt = -2*np.pi*np.cos(2*np.pi*(x-t))
  gx = 2*np.pi*np.cos(2*np.pi*(x-t))
  return gt + g1*gx

#パラメータの設定
N = 201
T = 0.5
mu = 1
num = 40

#非線形保存則の計算
x, tn, sol = conservationlaw_cf(flux1, initial1, righthand1, N, T, mu, num)

#図の表示
plot_solution(x, sol)

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

In [8]では，$N=301$としているが，近似の精度はあまり良くないようである． そこで，誤差の挙動を調べるためにIn [9]を用意する．

In [9]: 
#誤差の計算、非線形保存則
def error_conservationlaw(flux, initial, righthand, exact_sol, T, mu):
  L = 1
  N = 40
  num = 40

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

  for k in range(kmax):
    N = N*2
    x, tn, sol = conservationlaw_cf(flux, initial, righthand, N, T, mu, num)
    err = 0.0
    for n in range(tn.shape[0]):
      tval = tn[n]
      err = np.maximum(err, np.linalg.norm(sol[n,:] - exact_sol(x, tval), ord=np.inf))
    h = L/N
    hv[k] = h
    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 [10]を実行することで，CL差分スキームが収束していることが確かめられる．すなわち，近似の精度はあまり良くないのは，$h$ や $\tau$ が十分には小さくとれていないことが原因であることがわかる．

In [10]:  
def initial1(x):
  return np.sin(2*np.pi*x) + 1.2

def flux1(s):
  return 0.5*s*s

def righthand1(x,t):
  g0 = 1.2
  g1 = np.sin(2*np.pi*(x-t)) + g0
  gt = -2*np.pi*np.cos(2*np.pi*(x-t))
  gx = 2*np.pi*np.cos(2*np.pi*(x-t))
  return gt + g1*gx

def solution1(x, t):
  g0 = 1.2
  g1 = np.sin(2*np.pi*(x-t)) + g0
  return g1

#パラメータの設定
T = 0.5
mu = 1.0

#HJ方程式の計算
hv, ev, rate = error_conservationlaw(flux1, initial1, righthand1, solution1, T, mu)

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

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


Out [10]: 
0.00625000, 0.663
0.00312500, 0.712
0.00156250, 0.797
0.00078125, 0.816
0.00039063, 0.880

Hamilton--Jacobi方程式に対する差分スキーム

次の形のHamilton--Jacobi方程式を $\mathbb{T}=\mathbb{R}/\mathbb{Z}$ 上で考える： \begin{align} u_t+H(u_x)&=g(x,t),&& x\in\mathbb{T},\ t\in (0,T),\tag{C.4a} \\ u(x,0)&=u_0(x),&& x\in\mathbb{T}.\tag{C.4b} \end{align} ただし， $H:\mathbb{R}\to \mathbb{R}$, $g:\mathbb{T}\times [0,T]\to \mathbb{R}$, $u_0:\mathbb{T}\to \mathbb{R}$ は与えられた関数である．講義では，$g\equiv 0$の場合のみを説明したので注意せよ．

次のように記号を定める． \begin{equation*} \delta^{+}_h u_i^n=\frac{u^n_{i+1}-u_i^n}{h},\quad \delta^-_h u_i^n=\frac{u^n_{i}-u_{i-1}^n}{h}. \label{eq:fdq1} \end{equation*} このとき，(C.4a)に対するLax–Friedrichs (LF) 差分スキームは， \begin{align} u_i^{n+1}&=\frac{u_{i-1}^n+u_{i+1}^n}{2}-\tau H\left(\frac{\delta^-_h u_i^n+\delta^+_h u_i^n}{2}\right)+\tau g_i^n \\ &=u^n_i-\tau \left[H\left(\frac{\delta^-_h u_i^n+\delta^+_h u_i^n}{2}\right)-\frac{h}{\tau}\frac{\delta_h^+u_{i}^n-\delta_h^-u_i^n}{2}\right]+\tau g_i^n \end{align} である．

風上差分スキームを考える際には，次の形の特別な $H$ のみを扱う： \begin{equation} \exists \hat{p}\in\mathbb{R},\quad H'(\hat{p})=0, \quad \begin{cases} H'(p)\le 0 & (p\le \hat{p}),\\ H'(p)\ge 0 & (p>\hat{p}). \end{cases} \label{eq:HJ1d} \end{equation} このとき，HJ方程式に対する風上差分スキームは， \begin{equation} u_i^{n+1}=u_i^n-\tau \left[ \theta(\delta^+_h u_i^n)(H(\delta^+_h u_i^n)-H(\hat{p})) \right. \\ \left. +(1-\theta(\delta^-_h u_i^n))(H(\delta^-_h u_i^n)-H(\hat{p})) +H(\hat{p})\right]+\tau g_i^n \end{equation} である．ただし， \[ \theta (p)= \begin{cases} 1 & (s\le \hat{p})\\ 0 & (s> \hat{p}) \end{cases} \] としている．これらは，[CL1984]による（[FF2014]も参考になる）．

用語の確認をしておこう．$g\equiv 0$ とする．次の形で表現できる差分スキームを，差分形式の差分スキームと言う： \begin{equation} u_i^{n+1}=u_i^n-\tau \varphi(\delta^-_h u_i^n,\delta^{+}_h u_i^n). \end{equation} そして，$\varphi$ を数値Hamiltonianと言う． LF差分スキームと風上差分スキームは差分形式の差分スキームである．さらに，差分形式の差分スキームを \[ u_i^{n+1}=\Phi(u_{i-1}^n,u_i^n,u_{i+1}^n) \] の形に書いたとき，それが $[-R,R]$ で単調であるとは，$v_1,v_2,v_3\in\mathbb{R}$ が， $|\delta_h^- v_2|\le R$，$|\delta_h^+ v_2|\le R$ を満たす限りにおいて，$\Phi(v_1,v_2,v_3)$ が各 $v_j$ $(j=1,2,3)$ について非減少関数となるときを言う．単調な差分形式の差分スキームが，単調性，安定性，順序保存性などとても良い性質を持つことは講義で説明した．

重要なことのみ簡単に復習する（いずれも，$g\equiv 0$ の場合のみを考えている）：

1. $R>0$ とする．F差分スキームは，CFL条件 \begin{equation} \tag{C.5} \lambda \max_{|p|\le R}|H'(p)|\le 1 \end{equation} の下で，$[-R,R]$ で単調である．
2. $Q=\mathbb{T}\times [0,T]$ とおき，$R=\|u_x\|_{L^\infty(Q)}+1$ とする．ただし，$u$ は(C.4)の滑らかな解である． 初期値を $|\delta^+_hu_i^0|\le R$ となるようにとると，LF差分スキームの解 $u_i^n$ は，CFL条件 (C.5)の下で，$|\delta^+_hu_i^n|\le R$ と満たし， \[ \max_{(x_i,t_n)\in Q}|u(x_i,t_n)-u_i^n|\le Ch \] が成り立つ．ここで，$C$ は $u$，$H$，$R$，$T$ に依存して定まる正定数である（$h$ や $\tau$ には依存しない） ．
3. $R>0$ とする．風上差分スキームは，CFL条件 \begin{equation} \lambda \max_{|p|\le R}|H'(p)|\le \frac{1}{2} \end{equation} の下で，$[-R,R]$ で単調である．

LF差分スキームを実行するためには，CFL条件(C.5)に注意を払わなければならない．しかし，数値計算の際に$R=\|u_x\|_{L^\infty(Q)}+1$ の値を正確に取得することは難しい．そこで，次の様に考える．初期関数の導関数について $|\frac{d}{dx}u_0(x)|$ の最大値 $R_0$ を計算して，$R=R_0+1$ と設定する． その上で， $\Delta p= 2R/1000$，$p_j= -R+j\Delta p$ $(j=0,\ldots,1000)$ とし， \[ \left|\frac{H(p_{j+1})-H(p_j)}{\Delta p}\right|\qquad (j=0,\ldots,1000-1) \] の最大値 $\mathrm{d}H_{\max}$ を取得する．これを計算する際にはIn [5]を微修正したIn [11]を利用する．

In [11]:  
def find_maximum_dH(HH, dinit):
  p1 = np.max(np.abs(dinit))
  Np = 1000
  xp = np.linspace(-p1, p1, Np)
  yp = HH(xp)
  dHp = (yp[1:] - yp[:-1])*Np/(2*p1)
  return p1, np.max(np.abs(dHp))

はじめに，$\mu=\lambda \max_{|p|\le R}|H'(p)|$ の値を，$0 < \mu\le 1$ の範囲で指定しておく． そして，$h$ を $h=1/N$ で定めた後に，$\tau=\mu h /\mathrm{d}H_{\max}$ および $\lambda=\tau/h$ とする．それが，In [12]である．

In [12]:  
def hamiltonjacobi_lf(Hamil, initial, func_g, N, T, mu, num):
  # x-grids
  L = 1.0
  h = L/N
  xgrid = np.linspace(0.0, L, N)
  xtmp = xgrid[:-1].copy()

  # initial setteings
  tgrid = np.array([0.0])
  u = initial(xtmp)
  sol = np.copy(np.append(u,u[0]))
  
  # CFL condition
  uu = initial(xgrid)
  du = (uu[1:] - uu[:-1])/h
  pmax, abs_dH = find_maximum_dH(Hamil, du)

  # t-grids
  tau = mu*h/abs_dH
  lam = tau/h
  nmax = int(T/tau)
  step = int(max(1, nmax/num))

  # iteration 
  tnow = 0.0
  for n in range(nmax):
    coef_g = func_g(xtmp, tnow)
    u = 0.5*(np.roll(u,1) + np.roll(u,-1)) - tau*Hamil(0.5*(np.roll(u,-1)-np.roll(u,1))/h) + tau*coef_g
    tnow = (n+1)*tau
    if (n+1)%step==0:
      sol=np.vstack((sol,np.append(u,u[0])))
      tgrid=np.append(tgrid, tnow)

  return xgrid, tgrid, sol

ハミルトニアン，初期関数と右辺の関数を， \begin{align} H(p)&=\frac12 p^2,\\ u_0(x)&=\sin(2\pi x),\\ g(x,t)& =-2\pi \cos(2\pi(x-t))+2\pi^2\cos^2(2\pi(x-t)) \end{align} と選ぶ．すると，(C.4)の厳密解は， \begin{equation} u(x,t) =\sin(2\pi (x-t)) \end{equation} となる．これらに対して， In [13] のようにLF差分スキームを実行してみる． 結果を，図C.4とC.5に示す．$N=301$としているが，近似の精度はあまり良くないようである．誤差の挙動を調べることで，収束自体はしていることが確かめられる．すなわち，近似の精度はあまり良くないのは，$h$ や $\tau$ が十分には小さくとれていないことが原因であることがわかる．

In [13]:  
# Hamiltonian
def Hamil(p):
  return 0.5*p*p

# 初期値
def initial(x):
  return np.sin(2.0*np.pi*x)

#右辺関数
def func_g(x, t):
  y = np.pi*np.cos(2.0*np.pi*(x - t))
  return 2.0*y*(y - 1.0)

#パラメータの設定
N = 301
T = 1.0
mu = 1.0
num = 40

#HJ方程式の計算
x, tn, sol = hamiltonjacobi_lf(Hamil, initial, func_g, N, T, mu, num)

#図の表示
plot_solution(x, sol)

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

最後に誤差を観察しよう．そのために， In [14]を用意して In [15]のように実行する． 結果は，Out [15]のようになる．この結果からは，誤差は $h$ に比例して $0$ に減衰する．

In [14]:  
#誤差の計算、HJ方程式
def error_hamiltonjacobi(Hamil, func_g, initial, exact_sol, T, mu):
  L = 1.0
  N = 40
  num = 40

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

  for k in range(kmax):
    N = N*2
    x, tn, sol = hamiltonjacobi_lf(Hamil, initial, func_g, N, T, mu, num) 
    err = 0.0
    for n in range(tn.shape[0]):
      tval = tn[n]
      err = np.maximum(err, np.linalg.norm(sol[n,:] - exact_sol(x, tval), ord=np.inf))
    h = L/N
    hv[k] = h
    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]: 
# Hamiltonian
def Hamil(p):
  return 0.5*p*p

# 初期値
def initial(x):
  return np.sin(2.0*np.pi*x)

#右辺関数
def func_g(x, t):
  y = np.pi*np.cos(2.0*np.pi*(x - t))
  return 2.0*y*(y - 1.0)

def solution(x, t):
    return np.sin(2.0*np.pi*(x-t))

#パラメータの設定
T = 1.0
mu = 1.0

#HJ方程式の計算
hv, ev, rate = error_hamiltonjacobi(Hamil, func_g, initial, solution, T, mu)

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

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

Out [15]: 
0.00625000, 1.007
0.00312500, 0.998
0.00156250, 0.999
0.00078125, 0.992
0.00039063, 1.000

問題

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

課題

課題C.1

例題を自分で設定して，Hamilton--Jacobi方程式に対して，風上差分スキームを実行せよ．（計算結果を示すのみで良い．）

参考文献

• [CL1984] M. G. Crandall and P. L. Lions, Two approximations of solutions of Hamilton--Jacobi equations, Math. Comp. 43, No. 167 (1984) 1--19.
• [FF2014] M. Falcone and R. Ferretti, Semi-Lagrangian Approximation Schemes for Linear and Hamilton--Jacobi Equations, Society for Industrial and Applied Mathematics (SIAM), 2014.