$ \newcommand\od[3]{\frac{d^{#3} {#1}}{d {#2}^{#3}}} \newcommand\pd[3]{\frac{\partial^{#3}{#1}}{\partial {#2}^{#3}}} $
　

2021年9月4日 洛星高校 土曜講座 Part2

最終更新：2021/09/14

土曜講座の内容（微分方程式のシミュレーション）を公開します. 多くを板書で説明したので要約となっていますが, （特に数式については）できるだけ詳しく書くようにします.

正直高校生レベルじゃなくなりました. 雰囲気楽しんでもらえたらと思います.

その他の計算例

習うより慣れろで色んな方程式をみてみよう. 以下, $c$ は正の定数とします.

波動方程式

波動方程式は以下のように表されます.

\begin{align}\large{ \frac{\partial^\color{red}{2} u}{\partial t^\color{red}{2}} = c \pd{u}{x}{2} \label{waveeq} }\end{align}

熱方程式との違いは, 赤字で示したように時間微分が2階となっているだけです. これが解にどう影響するでしょうか?

早速離散化して計算してみましょう. Part1の熱方程式の離散化を参考にすると, 次のように離散化されることがわかります.

\begin{align}\large{ \frac{u_{i, j+1} - 2u_{i, j} + u_{i, j-1}}{{\Delta t}^2} = c \frac{u_{i+1, j} - 2u_{i, j} + u_{i-1, j}}{{\Delta x}^2}. }\end{align}

プログラミングしやすい形に変形して,

\begin{align}\large{ u_{i, j+1} = 2u_{i, j} - u_{i, j-1} + c \frac{{\Delta t}^2}{{\Delta x}^2} (u_{i+1, j} - 2u_{i, j} + u_{i-1, j}). }\end{align}

とします. また, 熱方程式のときのように初期状態を決める必要があるので, 以下のようにしましょう. これはただのサイン関数です. また, 境界条件として両端は0に固定します.

name

どんな結果になると思いますか? 以下に結果を示します.


上図のように固定端で振動するひも, 縄跳びのようになります. 波動方程式という名前も納得です. なぜこのような結果になるのでしょうか?

参考として動画の下には $u$ の空間方向の二階微分の値を示しています. 熱方程式の場合は, この値は $u$ の時間に関する1回微分, 則ち変化速度に比例するのでした. 一方, 波動方程式では $u$ の時間に関する2回微分, 則ち加速度に比例します. なので, グラフの形が上に凸であるかぎり, $u$ は減少方向に加速され, 0 に達しても熱方程式のようにとまることなく, 最高速度で通り越してしまうのです.

コードをいかにおいておきます.

MATLABコード（クリック） MATLAB 
L    = 1;
DT   = 0.05;
DX   = 0.05;
ITER = 200;
C    = 0.7^2;

% preprocessing
R2 = C * DT^2/DX^2;
X  = 0:DX:L;
U0 = sin(pi*X)/5;
NX = length(0:DX:L);
u  = zeros(ITER, NX);
u  = [U0; U0; u];

% main
for j = 2:ITER
    for i = 2:NX-1
        u(j+1, i) = R2*u(j, i+1) + 2*(1-R2)*u(j, i) + R2*u(j, i-1) - u(j-1, i);  
    end
end

移流方程式

次に, 以下のような方程式をみてみましょう. 熱方程式と比べると, 両辺一回微分で, 簡単になりました.

\begin{align}\large{ \pd{u}{t}{} = -c \pd{u}{x}{} \label{advection} }\end{align}

とりあえず離散化してみよう.

\begin{align}\large{ \frac{u_{i, j+1} - u_{i, j}}{{\Delta t}} = -c \frac{u_{i, j} - u_{i-1, j}}{{\Delta x}}. }\end{align}

プログラミングしやすい形に変形して,

\begin{align}\large{ u_{i, j+1} = u_{i, j} - c \frac{{\Delta t}}{{\Delta x}} (u_{i, j} - u_{i-1, j}) }\end{align}

を計算します. 初期条件は以下のようにしましょう.

name

結果は以下のようになります.


このように, どのような初期状態を与えても, その状態が速さ $c$ で流れていくような結果となります. これが, 移流方程式と呼ばれる所以です. なお, 高さがちょっとずつ低くなっているのは計算による誤差なので, 無視して下さい. この点については, 次の節で少し触れます.

さて, なぜこのような結果になるのか, ですが, 下段のグラフを見て下さい. これは移流方程式の右辺の値（ただし, $c$ は定数なので無視）です. つまり, この値に比例して $u$ の値は時間的に変化します.

変化量が正の値で大きい（つまり傾きが多きい）ところでは, 大きな速さで値が減少します. 逆に, 変化量が小さい（傾かが小さい）と, 小さな速さで変化します. その結果, あたかも波形が保たれたまま一方向の流れていくようにみえるのです.

コードを以下においておきます.

MATLABコード（クリック） MATLAB 
L    = 1;
DT   = 0.01;
DX   = 0.01;
ITER = 600;
C    = 0.2;

% preprocessing
R = C * DT/DX;
X   = 0:DX:L;
U0  = exp(-200*(X-0.15).^2)/5;
NX = length(0:DX:L);
u = zeros(ITER, NX);
u = [U0; u];

% main
for j = 1:ITER
    for i = 2:NX
        u(j+1, i) = u(j, i) - R*(u(j, i) - u(j, i-1)); 
    end
end

Burgers方程式

このあたりからかなり発展的な内容です. （非粘性）Burgers方程式と呼ばれる方程式は以下のように記されます.

\begin{align}\large{ \pd{u}{t}{} = -\color{red}{u} \pd{u}{x}{} \label{burgers} }\end{align}

移流方程式のときは, $\color{red}{u}$ の部分が 正の定数 $c$ でした. しかし, こんどは $u$ 自身が入っており定数ではありません. この小さな変化がどんな結果を生むでしょうか? とりあえず離散化してみよう.

\begin{align}\large{ \frac{u_{i, j+1} - u_{i, j}}{{\Delta t}} = -u_{i-1, j} \frac{u_{i, j} - u_{i-1, j}}{{\Delta x}}. }\end{align}

プログラミングしやすい形に変形して,

\begin{align}\large{ u_{i, j+1} = u_{i, j} - u_{i-1, j} \frac{{\Delta t}}{{\Delta x}} (u_{i, j} - u_{i-1, j}) }\end{align}

を計算します. 初期条件は移流方程式のときと同じにします.

結果は以下のようになります.


このへんになると, 僕の専門とはいよいよ離れて来るのであまり詳細を語れませんが, このようになる理由を簡単に説明します.

先程の移流方程式では, $\color{red}{u}$ の部分が 正の定数 $c$ で, 速さ $c$ で波が移動していくのでした. それが, 今回は $u$ となっており, これは $u$ の値が大きい（即ちグラフの高さが高い）ほど波は速く移動し, 逆に小さいと移動速度も小さいということです. なので, 今回のような初期状態を与えると, 山のてっぺんはどんどん右に移動しようとするものの, 裾野はゆっくり進むので, 波は前に突っ立ったような形になります. 実際, このような現象を突っ立ちという.

現実現象としては, 海で突っ立ったあとに波が崩れるようなものが馴染み深いと思う（この方程式で表現しきれるわけではないけれど）. これは砕波とよばれるもので, ぜひ画像検索などしてみてほしい.

コードを以下においておきます.

MATLABコード（クリック） MATLAB 
L    = 1;
DT   = 0.01;
DX   = 0.01;
ITER = 600;
C    = 0.2;

% preprocessing
R = C * DT/DX;
X   = 0:DX:L;
U0  = exp(-200*(X-0.15).^2)/5;
NX = length(0:DX:L);
u = zeros(ITER, NX);
u = [U0; u];

% main
for j = 1:ITER
    for i = 2:NX
        u(j+1, i) = u(j, i) - u(j, i-1)*(u(j, i) - u(j, i-1)); % 移流方程式と比べて変わったのはここだけ!
    end
end

KdV方程式

これは, 完全に発展レベルを超えていますが, 見て楽しめると思うので紹介します. 次の方程式を, Korteweg–De Vries（KdV）方程式といいます. 移流方程式と異なる点は, $u$ の $x$ に関する三階微分の項が追加されていることです. これがどう働くでしょうか?

\begin{align}\large{ \pd{u}{t}{} = -6u \pd{u}{x}{} - \pd{u}{x}{3} \label{kdv} }\end{align}

これは, 普通に離散化してもあまりに効率が悪いので, フーリエ変換を使った方法で計算しました. もし興味ある人は, Numerical Solution of the KdVなどを見て下さい.

結果は以下のようになります.


不思議な動きですね. 波の突っ立ちが起きそうになった瞬間波が分裂してしまいました. なぜこのようになるのかを私も完璧に理解している訳ではありませんが, 参考として下段にKdV方程式の右辺第一項, 第二項の時間変化も示しています.

見てもらいたいのは最下段, 三階微分です. 三階微分の値は波が滑らかな間は小さな値ですが, 波が突っ立つにつれて急激に大きくなり, 一回微分項を打ち消すように値をとります. これが, 波の突っ立ちを防ぎ, また, 波を孤立させるのです. こういった波を孤立波, ソリトンといいます. そして, 孤立した波の背部でまた突っ立ちが起きそうになるものの, それを防ぐように三階微分の値が大きくなり, 波の孤立の連鎖が起きています.

この記事では, これ以上深く突っ込みませんが, 物理・数学的にはかなり奥が深い内容であるようです. 興味がある人は, KdV方程式, ソリトンなどで調べてみて下さい.

最後に, もう少し先まで計算した結果を以下に載せておきます.

数値計算の難しい側面

さて, ここまで色々な偏微分方程式をみて, 数式とシミュレーション結果の関係を説明してきました. 少しでも数式の意味が想像出来るようになって貰えれば幸いです. また, 離散化することで簡単にシミュレーション出来ることも示しました.

しかし, 実際シミュレーションはそれほど簡単にうまくいく訳ではありません. むしろ上手くいかないのが普通で, 今回の記事もうまくいく計算条件を探して紹介したのが実情です. こんな単純な離散化でどんな数式も計算できるなら第3の科学と呼ばれることもないでしょう. ここでは, 少し熱方程式を取り上げて, 数値計算の難しい側面をお見せしようと思います.

ここでは, 発散と実用的な計算時間について述べます.

解の発散

離散化して計算する際, $\Delta t$ や $\Delta x$ を決める必要がありました. それぞれの値を大きくとると, 粗い計算となり, 小さくなると高い精度の計算となることが想像できるとおもいます. 例えば, 熱方程式の例では, $\Delta t = 0.002$, $\Delta x = 0.2$ として計算しました. 今, もう少し空間的に高精度で計算したいとしましょう. その場合, $\Delta x = 0.2$ を小さくすればよく, 例えば $\Delta x = 0.15$ としてみましょう. すると結果は以下のようになります.


値が振動し, 極めて大きな値となっていることがわかります. このような計算の破綻を発散といいます. 実は, 計算の解像度を上げたいからといって $\Delta x$ を無闇に小さくしても良い計算結果が得られるとは限らないのです. 理由は複雑なので割愛しますが, 興味がある人はCFL条件などで調べてみて下さい.

実用的な計算時間

実は, 今回の場合は $\Delta x$ を小さくするなら, $\Delta t$ も十分に小さくすれば計算がうまくいきます. しかし, $\Delta x$ を1/2倍にするときに, $\Delta t$ も1/2倍にすれば良いという単純な話ではありません. 今回の場合, 理論的に $\Delta t$ は 1/2.8 倍にする必要があります. これが計算時間にどう影響するでしょうか?

まず, $\Delta x$ を1/2倍にするときに, 2倍の計算量が必要であることを確認しておきます. なぜなら, これは2cm刻みで計算していたシミュレーションを, 1cm刻みで計算するようなことに相当するからです. 時間刻みを小さくするときも同様です. よって, $\Delta x$ を1/2倍にするときには, 最初の精度で計算時間が1分だとすると 1[分]×2×2.8で5.6分の計算時間がかかります.

もっと高精度にする場合をみてみましょう. たとえば, $\Delta x$ を1/10倍にすると $\Delta t$ は1/70倍以上にする必要があります. この場合, 1×10×70で700分以上かかります. このように, 空間解像度を $c$ 倍すると, 計算時間は $c$ 倍以上, $c$ 倍をはるかに超える時間がかかるのが一般的です. 二次元, 三次元の計算になるとこの増分はもっと大きくなります.

つまり, 我々の計算に費やせる時間が有限であることを考えると, 必ずしも自分が必要な精度でシミュレーションをすることができるわけではないのです. こういったときには, より効率のよい計算方法, 離散化方法を考える必要があり, 複雑な数学が必要となってきます. このような数学の分野は, 数値解析などと呼ばれております. 実際, 上の例でKdV方程式を計算する際は, 単純な離散化を諦めてフーリエ変換を用いた方法を使わせて頂いています.

この分野は私の専門ではないので, 多くは語れませんが, 東京大学松尾教授ウェブサイトの紹介文が非常にわかり易く, 魅力的です. ぜひ一読してみて下さい.

その他

他にも, 面白い例をわかり易く書けそうでしたら, ここに今後も追加していこうと思います.

この記事のTOP    BACK    TOP