数値解析の基礎：10分解説

数値解析のアプローチ

• 自然現象や社会現象など世の中の多くの現象は微分方程式によって表すことができる（当然その限りではない）。
• 数値解析とは微分方程式によってモデル化された現象をコンピュータで取り扱える形にして計算を行なうことを示す。
• 微分方程式を離散化することにより代数方程式を導出する。離散化手法のうち
  、流体解析でよく使われる差分法について説明する

※ 茨城大学 田中伸厚先生の講義資料を引用 : http://www.mech.ibaraki.ac.jp/~tanaka/Lecture/sim1/part1.pdf

離散化について -差分法を例に-

ある現象が次の常微分方程式で表すことができるとする。

:定数

上式は連続系でありそのままコンピュータで取り扱うことができない。

そこで上式を

のようにコンピュータ上で取り扱える式（代数方程式＝四則演算のみの式）に変換することを離散化（差分法）という。

また、代数方程式は添字を用いて以下のように記述することが多い。以降は以下の表記を使用する。

計算例1：放射性物質の崩壊

現象

放射性物質の原子核数の時間変化

出典：https://studyphys.com/half-life/

数学モデル（微分方程式）

以下の常微分方程式で表すことができる。

:崩壊定数

数値モデル（代数方程式）

差分法を用いると以下のとおりとなる。

この式よりの値（初期条件）を与えて逐次計算を行うことにより数値解を得ることができる。

厳密解

前出の常微分方程式は厳密解（数値計算を用いずに解析的に得られる解。解析解とも呼ばれる）を求めることができる。

導入過程は省略するとが、変数分離を用いて式変形することにより下式が得られる。

本式による計算結果と数値計算の結果を比較することにより計算の妥当性を評価することができる。

計算結果

導出した代数方程式を用いて数値解析を実施する。
を0.1,0.5,1.0,1.5と変化させた場合の計算結果と厳密解との比較を示す。
計算条件は、は100、は1とする。

import numpy as np
import matplotlib.pyplot as plt

N0 = 100
k = 1.0
exact = lambda T: N0 * np.exp(-k*T)

def simulation(T):
    N = np.zeros_like(T)
    N[0] = N0
    for i in range(1, len(N)):
        N[i] = N[i-1] - dt*k*N[i-1] 
    return N

def mkfig(ax, T):
    ax.plot(T, exact(T),label='exact')
    ax.plot(T, simulation(T),label='simulation')
    ax.set_xlabel('T')
    ax.set_ylabel('N')
    ax.legend(loc = 'upper right') 

の計算結果

dt = 0.1
T = np.arange(0, 20, dt)
fig, ax = plt.subplots()
mkfig(ax, T)

の計算結果

dt = 0.5
T = np.arange(0, 20, dt)
fig, ax = plt.subplots()
mkfig(ax, T)

の計算結果

dt = 1.0
T = np.arange(0, 20, dt)
fig, ax = plt.subplots()
mkfig(ax, T)

の計算結果

dt = 1.5
T = np.arange(0, 20, dt)
fig, ax = plt.subplots()
mkfig(ax, T)

• によって解析結果が異なる。
• は、厳密解と数値解析による差はほとんどないが、が大きくなるにつれて差が大きくなる。は異常な結果を示している。
• 差分法ではを小さく設定することが望ましいが、計算資源の制約で限界がある。

他の差分法

• 今回はもっと基本的な前進差分スキーム（nの値からからn+1の値を求める）を使用したが、後退差分スキーム（n+1の値からからn+1の値を求める）も存在する。
• 時間発展の前進差分スキームは陽解法、後退差分スキームは陰解法と呼ばれる。

• 以下に示す計算結果のとおり、でもそれなりの数値計算結果が得られることがわかる。
• 陰解法はメリットしか無いように思われるが、陰解法では上記のような単純な式形になることは稀であり、ほとんどのケースで複雑な式形になり、数値解を得るために収束計算が必要となる。

import numpy as np
import matplotlib.pyplot as plt

N0 = 100
k = 1.0
exact = lambda T: N0 * np.exp(-k*T)

def simulation2(T):
    N = np.zeros_like(T)
    N[0] = N0
    for i in range(1, len(N)):
        N[i] = N[i-1]/(1+dt*k)
    return N

def mkfig2(ax, T):
    ax.plot(T, exact(T),label='exact')
    ax.plot(T, simulation2(T),label='simulation')
    ax.set_xlabel('T')
    ax.set_ylabel('N')
    ax.legend(loc = 'upper right') 

の計算結果

dt = 1.5
T = np.arange(0, 20, dt)
fig, ax = plt.subplots()
mkfig2(ax, T)

計算例2：移流方程式

現象

移流とは下図のようにある物理量が流れによって運ばれる現象を示す。この現象を計算対象とする。

数学モデル（微分方程式）

以下の偏微分方程式で表すことができる。

:速度

数値モデル（代数方程式）

今回は3つの差分スキームを使用する。

前進差分

後退差分

中心差分

この式よりの値（初期条件,空間分布を持つ）を与えて逐次計算を行うことにより数値解を得ることができる。
なお、境界条件も必要であるが今回は境界付近を計算しないため0とする。

厳密解

一次元の移流方程式は条件によって厳密解を得ることができる。
を与えると厳密解として次式が得られる。

計算結果

導出した代数方程式を用いて数値解析を実施する。
計算条件は、、、として、厳密解との比較を行なう。

import numpy as np
import matplotlib.pyplot as plt

c = 1.0
X = np.linspace(-5,5,101)
dx = 0.1
exact = lambda t: np.exp(-(X-c*t)**2)

def s01(dt): #前進差分
    unew = exact(0)
    for _ in range(int(2/dt)):
        u = unew.copy()
        for i in range(1,len(u)-1):
            unew[i] = u[i] - c * dt/dx * (u[i+1]-u[i])
    return unew

def s02(dt): #後退差分
    unew = exact(0)
    for _ in range(int(2/dt)):
        u = unew.copy()
        for i in range(1,len(u)-1):
            unew[i] = u[i] - c * dt/dx * (u[i]-u[i-1])
    return unew

def s03(dt): #中心差分
    unew = exact(0)
    for _ in range(int(2/dt)):
        u = unew.copy()
        for i in range(1,len(u)-1):
            unew[i] = u[i] - c * dt/dx * (u[i+1]-u[i-1])/2
    return unew

前進差分

適切な解を得ることができない。

後退差分

波形が拡散してピークが減衰する。

中心差分

若干の差異はあるが他の差分スキームと比較すると厳密解に近い結果になっている。

• 差分スキームによって解析結果が異なる。
• 前進差分はをいくら小さくしても適切な解を得ることができない。
• 後退差分は拡散するため、を十分に小さくとる必要がある。
• 今回のケースでは中心差分が良さそうであるが、適切な解が得られない場合が多い。（発散し易い。）
• その場合、より高精度なスキームが必要となり代数方程式の導出が複雑となる。⇒ 不定流計算などでは煩雑な代数方程式が必要

陰解法

一次元移流方程式の陰解法は前出の常微分方程式と比べてかなり複雑な式形となる。

例として陰解法で中心差分の代数方程式を導出すると次式となる。

上式のみでは、の項の値を求めることができない。そのため、他の空間離散点（など）の式を全て連立させた次式より解を計算する。

左辺の対角行列の逆行列を計算することにより、の項を求めることができる。

ここでは計算は省力するが、陰解法は陽解法と比べて煩雑な計算が必要となることが理解できる。

補足として、上式の左辺の1つ目の行列は三重対角行列となっており、TDMA法を用いて少ない計算ステップ数で逆行列を求めることができる。この方法は不等流計算の平均流速公式レベル3の解法でも使用する。

離散化方法のポイント

一般的に精度と安定性は両立しない。目的に応じて差分スキームを選択する必要がある。

精度

• 離散化による打ち切り誤差の大きさ(○次精度のスキームなどと呼ぶ)

安定性

• 計算誤差が時間の経過とともに発展しないこと。
  • 中心差分は計算誤差が発展しやすく不安定なスキームである。
  • 陰解法は安定性が高い。

まとめ

• 数値解析は微分方程式を差分法などにより離散化した代数方程式から数値解を得ることである。
• 差分法では、の大きさによって解析結果が大きく異なるため、適切な条件設定が必要である。
• 差分法には様々な方法（スキーム）があり、目的に応じた選択が必要である。

⇒ 今回の勉強で実施する不等流計算、貯留関数法では、差分スキームはそれほど重要ではないがの大きさについてはある程度考慮する必要がある。

⇒ 不定流計算などを実施する場合は、差分スキームが重要