Pythonで2次元ラプラス方程式を数値計算する方法

  • このエントリーをはてなブックマークに追加

時間によって変化しない定常状態を示すラプラス方程式は様々な物理現象の基礎となっています。ここではラプラス方程式の差分化から説明します。また、学習の理解を深めるために簡単な2次元のラプラス方程式をPythonを使って数値計算する方法を紹介します。

こんにちは。wat(@watlablog)です。ここでは2次元のラプラス方程式を離散化して数値計算する方法を紹介します

ラプラス方程式の概要

基礎方程式

ラプラス方程式Laplace's Equation)とは、式(1)に表す楕円型偏微分方程式です。

\[ \nabla ^{2} \phi = 0 \tag{1} \]
Advertisements

ここで\(\phi\)は任意の物理量、\(\nabla ^{2}\)はラプラス作用素Laplace operator)またはラプラシアンLaplacian)と呼ばれ\(^{[1]}\)、空間2階微分を意味します。2次元空間\(\mathrm{R}^{2}\)であれば式(1)のラプラス方程式は式(2)と書くことができます。

\[ \frac{\partial^{2} \phi}{\partial x^{2}} + \frac{\partial^{2} \phi}{\partial y^{2}} = 0 \tag{2} \]

この方程式は時間項が含まれていなく「定常状態または熱平衡状態」を表現しており、例えば境界に温度や電圧といったスカラー量が与えられた時、十分に時間が経過した後のスカラー分布を意味しています。

ラプラス方程式は電磁気学、流体力学、天文学といった様々な分野で多用されています。

今回筆者は計算力学技術者試験(熱流体)の学習でラプラス方程式を学んでいます。数値流体力学の分野ではナビエ・ストークス方程式圧力項を解く時にこの方程式を扱います(実際にはラプラス方程式の右辺が0でないポアソン方程式を解きますが)。

ここではポアソン方程式を学ぶ前に、より簡単なラプラス方程式の数値解析方法を確認します。

離散化

テイラー展開を使って2次精度中心差分式を作る

この記事では式(2)の2次元ラプラス方程式を例に数値計算を行います。
数値計算は既に以下の記事で扱ったように、離散化が必要です。

Pythonで1次元移流方程式を離散化して数値流体力学に入門する
Pythonで拡散方程式を数値計算してアニメーションを作成する

式(2)は空間項のみで構成されており、先の拡散方程式と同様に2次精度中心差分で離散化しましょう。

このような2変数の方程式でもやり方は同様です。まずは式(2)の中で\(x\)から考えます。

中心差分は前進方向後退方向を求めれば良さそうです。

2変数で表す任意の物理量を\(q(x, y)\)として、\(x\)方向に\(\Delta x\)だけずらした位置の前進方向の物理量\(q(x+\Delta x, y)\)はテイラー展開すると式(3)と求めることができます。

\[ q(x+\Delta x, y) = q(x, y) + \Delta x \frac{\partial q}{\partial x} + \frac{1}{2!}(\Delta x)^{2} \frac{\partial^{2} q}{\partial x^{2}} + \frac{1}{3!}(\Delta x)^{3} \frac{\partial^{3} q}{\partial x^{3}} + \cdots \tag{3} \]

同様に、後退方向の物理量\(q(x-\Delta x, y)\)は式(4)となります。式(3)と比べ奇数階微分項の符号が反転しています。

\[ q(x-\Delta x, y) = q(x, y) - \Delta x \frac{\partial q}{\partial x} + \frac{1}{2!}(\Delta x)^{2} \frac{\partial^{2} q}{\partial x^{2}} - \frac{1}{3!}(\Delta x)^{3} \frac{\partial^{3} q}{\partial x^{3}} + \cdots \tag{4} \]

式(3)と式(4)を加算すると、式(5)を得ます。ここで\(O((\Delta x)^{2})\)としたのは、2次以降の項を打ち切ったことを意味します。

\[ q(x+\Delta x, y) + q(x-\Delta x, y) = 2 q(x, y) + (\Delta x)^{2} \frac{\partial^{2} q}{\partial x^{2}} + O((\Delta x)^{2}) \tag{5} \]

式(5)の2階微分項を左辺、それ以外を右辺に移項し誤差項を省略して体裁を整えると、式(6)を得ます。この式が\(x\)方向の2次精度中心差分方程式です。

\[ \frac{\partial^{2} q}{\partial x^{2}} \approx \frac{q(x+\Delta x, y) - 2 q(x, y) + q(x-\Delta x, y)}{(\Delta x)^{2}} \tag{6} \]

\(y\)方向も同様に計算すると式(7)を得ます。

\[ \frac{\partial^{2} q}{\partial y^{2}} \approx \frac{q(x, y+\Delta y) - 2 q(x, y) + q(x, y-\Delta y)}{(\Delta y)^{2}} \tag{7} \]

ここまでの考え方で得た式(6)と式(7)を式(2)に代入すると、求めたい2次元ラプラス方程式の離散化形式を式(8)と得ることができます。

\[ \frac{\partial^{2} \phi}{\partial x^{2}} + \frac{\partial^{2} \phi}{\partial y^{2}} \approx \frac{q(x+\Delta x, y) - 2 q(x, y) + q(x-\Delta x, y)}{(\Delta x)^{2}} + \frac{q(x, y+\Delta y) - 2 q(x, y) + q(x, y-\Delta y)}{(\Delta y)^{2}} =0 \tag{8} \]

格子で考える

式(8)でも良いのですが、数値流体力学の一般的な格子点の座標表記に変換します。格子のイメージは下図です。\(i\)はプログラミングでよく使うので\(j, k\)としてみました。

格子の図

式(9)に格子点座標\(j, k\)で表現したラプラス方程式の2次精度中心差分方程式を示します。

\[ \frac{q_{j+1,k} - 2 q_{j,k} + q_{j-1,k}}{(\Delta x)^{2}} + \frac{q_{j,k+1} - 2 q_{j,k} + q_{j,k-1}}{(\Delta y)^{2}} =0 \tag{9} \]

漸化式

これまでの数値計算と同様に、漸化式を作って反復計算によりラプラス方程式の解を求めます。式(9)の\(q_{j, k}\)を左辺、それ以外を右辺にまとめることで式(10)の漸化式を得ます(演算自体は簡単ですが、結構面倒な計算かも)。

\[ \begin{eqnarray} \frac{q_{j+1,k} - 2 q_{j,k} + q_{j-1,k}}{(\Delta x)^{2}} + \frac{q_{j,k+1} - 2 q_{j,k} + q_{j,k-1}}{(\Delta y)^{2}} &=&0  \\ \frac{q_{j+1,k}}{(\Delta x^{2})} - \frac{2q_{j,k}}{(\Delta x)^{2}} + \frac{q_{j-1,k}}{(\Delta x)^2} + \frac{q_{j,k+1}}{(\Delta y)^{2}} - \frac{2q_{j,k}}{(\Delta y)^{2}} + \frac{q_{j,k-1}}{(\Delta y)^2} &=&0 \\ \frac{2q_{j,k}}{(\Delta x)^{2}} + \frac{2q_{j,k}}{(\Delta y)^{2}} &=& \frac{q_{j+1,k}}{(\Delta x^{2})} + \frac{q_{j-1,k}}{(\Delta x)^2} + \frac{q_{j,k+1}}{(\Delta y)^{2}} + \frac{q_{j,k-1}}{(\Delta y)^2} \\ \frac{2q_{j,k}((\Delta x)^{2}+(\Delta y)^{2})}{(\Delta x)^{2}(\Delta y)^{2}} &=& \frac{q_{j+1,k}}{(\Delta x^{2})} + \frac{q_{j-1,k}}{(\Delta x)^2} + \frac{q_{j,k+1}}{(\Delta y)^{2}} + \frac{q_{j,k-1}}{(\Delta y)^2} \\ q_{j,k} &=& \frac{(\Delta y)^{2}(q_{j+1,k}+q_{j-1,k})+(\Delta x)^{2}(q_{j,k+1}+q_{j, k-1})} {2((\Delta x)^{2}+(\Delta y)^{2})} \tag{10} \end{eqnarray} \]

ここで、式(10)は全て同時間ステップの値のみを使います。このようにラプラス方程式の漸化式は同一時刻における格子点の情報のみを使って反復計算を行うことで、分布を収束させる式となりました。

反復法による解法

式(10)の漸化式は同一時刻の情報のみを使って値を収束させていく形式となりましたが、式(10)は格子点の数(\((M-1) \times (N-1)\))分の連立方程式ができあがります。

多元連立方程式を反復法で解く方法として、当ブログで扱ったヤコビ法ガウス・ザイデル法SOR法が使えます。このページではその比較もしてみましょう。

それぞれの手法は以下の記事にまとめています。

Pythonで連立方程式をヤコビ法(反復法)で解く方法
Pythonで連立方程式をガウス・ザイデル法(反復法)で解く方法
連立方程式をSOR法で解くPythonコードと緩和係数のパラスタ

ここではそれぞれの手法における2次元ラプラス方程式の漸化式を示します。

ヤコビ法

式(11)に2次元ラプラス方程式をヤコビ法で解く場合の漸化式を示します。
添字\(n\)はステップ数を意味しています(時間ステップではありません)が、次のステップ\(n+1\)を更新するのに現在のステップ\(n\)のみを考慮するのが特徴です。また、 ヤコビ法は並列計算ができる手法でもあります。

\[ q^{n+1}_{j,k} = \frac{(\Delta y)^{2}(q^{n}_{j+1,k}+q^{n}_{j-1,k})+(\Delta x)^{2}(q^{n}_{j,k+1}+q^{n}_{j, k-1})} {2((\Delta x)^{2}+(\Delta y)^{2})} \tag{11} \]

ガウス・ザイデル法

式(12)に2次元ラプラス方程式をガウス・ザイデル法で解く場合の漸化式を示します。

ヤコビ法の右辺が全て\(n\)ステップ目の情報を使って漸化式を更新しているのに対し、ガウス・ザイデル法は既に求まった値(ここでは後退差分部分)を最新の値である\(n+1\)ステップ目の値を使って効率よく計算をしている特徴があります。

ガウス・ザイデル法はヤコビ法より少ない回数で収束に向かいますが、並列計算ができないというデメリットを持ちます。

\[ q^{n+1}_{j,k} = \frac{(\Delta y)^{2}(q^{n}_{j+1,k}+q^{n+1}_{j-1,k})+(\Delta x)^{2}(q^{n}_{j,k+1}+q^{n+1}_{j, k-1})} {2((\Delta x)^{2}+(\Delta y)^{2})} \tag{12} \]

SOR法

式(13)に示すSORSuccessive Over-Relaxation Method)は緩和係数\(\omega\)を使ってガウス・ザイデル法よりもさらに少ない回数で収束させようという手法です。基本的にはガウス・ザイデル法と同じ考え方のため、並列計算はできません。

\[ q^{n+1}_{j,k} = (1-\omega)q^{n}_{j,k} + \omega \frac{(\Delta y)^{2}(q^{n}_{j+1,k}+q^{n+1}_{j-1,k})+(\Delta x)^{2}(q^{n}_{j,k+1}+q^{n+1}_{j, k-1})} {2((\Delta x)^{2}+(\Delta y)^{2})} \tag{13} \]

SOR法の記事では、ヤコビ法の反復行列のスペクトル半径を行列から求めて最適な緩和係数\(\omega_{opt}\)を算出していましたが、格子サイズ\((m, n)\)がわかっている場合は式(14)でも求められるようです\(^{[2]}\)。

\[ \begin{eqnarray} \rho = \cos \frac{\pi}{m} + \cos \frac{\pi}{n}\\ \omega_{opt} = \frac{2}{1+\sqrt{1-(\frac{\rho}{2})^{2}}} \tag{14} \end{eqnarray} \]

初期場と境界条件

ラプラス方程式は任意の計算領域ドメイン)の中を反復計算します。そのため計算領域の端である境界には境界条件Boundary condition)が必要です。また、計算開始時には初期場(Initial field)として初期化(Initialization)も必要です。

今回は1辺が1の正方行列を100×100格子で用意し、上辺に境界条件1それ以外の3辺に0を与えました。そして境界以外の格子の値を0で初期化した条件を例にラプラス方程式を数値的に解いていきます。

以下はmatplotlibによるプロットの例です。ちなみに配列の表現はNumPyのスライス記法です。

境界条件と初期場

Pythonで2次元ラプラス方程式を解くコード

動作環境

このページでは以下の環境で動作確認を行ったコードを掲載しています。

Python環境

                
Python Python 3.9.6
PyCharm (IDE) PyCharm CE 2020.1
Numpy 1.21.1
matplotlib 3.4.3

PC環境

Mac OS macOS Catalina 10.15.7
CPU 1.4[GHz]
メモリ 8[GB]

ヤコビ法のコード

それではいよいよ2次元ラプラス方程式を数値計算で解くPythonコードを書いてみましょう。まずはヤコビ法の場合の全コードを以下に示します。

コードの構成としては、通常のシミュレーションソフトウェアをイメージして「初期場の設定」「境界条件の設定」「ソルバー」「ポスト処理結果の可視化)」という要素に分けて関数化しています。

コメント部分を順番に追って確認すればどのように計算しているかがわかると思います。今回は学習目的であり、式をそのまま書けるわかりやすさからダブルforループを書いていますのでプログラム実行速度は遅いです。

下図が結果です。よく見るラプラス方程式の結果がコンター図で確認できました。

100×100格子のラプラス方程式の数値解

プログラムは実行中に反復回数であるIterationをカウントしています。このプログラムではなんと2704回の反復回数を約34秒もかけて実行していました。反復回数はヤコビ法のアルゴリズム起因ですが、かかっている実行時間は単純にコーディングしたPythonのforループの遅さが目立っている結果と思います。

ガウス・ザイデル法のコード

続いて式(12)のガウス・ザイデル法のコードを以下に示します。ソルバー部分のみを書いています。古い分布を保持していますが計算のためではなく、収束計算のためです。その他は全て先ほどのコードと同様です。

結果は以下です。1930回の反復回数で約25秒でした。約1.4倍の速度UPができたようです(それでも遅いけど)。

SOR法のコード

最後は式(13)のSOR法で数値計算してみます。

先ほどと同様にソルバー部分のみのコードを示します。式(14)の最適な緩和係数を使って解くとどれだけはやくなるのでしょうか。

目を疑うほどの結果を得ました!SOR法はヤコビ法と比べ反復回数で23倍計算速度で18倍ほどの速度向上効果を得ました。

Pythonでforループを使っていますが体感時間的にもストレスを感じません。もっと大きい行列でもノートPCによる計算ができそうです。

ヤコビ法で2704回、ガウス・ザイデル法で1930回かかる計算がSOR法で117回とは…!おそるべしアルゴリズムの効果!

おまけ:収束仮定をアニメーションで確認してみる

ヤコビ法の収束過程動画

せっかくなので、当ブログ恒例の動画作成による可視化を行ってみましょう。
(コードは最後)

こちらはヤコビ法による反復計算中の解分布を可視化した動画です。ゆっくりと境界条件の値が伝播しているのがよくわかりますね。

ヤコビ法のアニメーション

ガウス・ザイデル法の収束過程動画

ガウス・ザイデル法の場合はこちら。動画のステップ数を先ほどより多めの20刻みにしていますが、基本的にはヤコビ法と同じような解変化をしていることがわかりました。

ガウス・ザイデル法のアニメーション

SOR法の収束過程動画

続いてSOR法の場合の動画です。こちらは何か間違ったか、と思うほど解の変化傾向が異なりました。端っこからどんどん最終形態に近づいていっている様子が面白いです。

SOR法のアニメーション

SOR法の場合のみ動画作成を含む全コードを紹介しておきます。是非お好みの境界条件を作ったりして遊んでみてください(境界条件は辺以外にもつくることができますよ!)。

まとめ

この記事は数値流体力学の中で、空間項のみで表されたラプラス方程式の概要を紹介しました。

ラプラス方程式は境界条件が決まると解が求まる楕円型偏微分方程式であることや、2次精度中心差分による離散化手法を説明しました。

式は一つ一つ丁寧に展開していったので、初学者にもやさしい内容になったのではと思っています(←)。

また、ここで紹介したラプラス方程式は連立方程式を反復法で解きましたが、当ブログで既に紹介したヤコビ法ガウス・ザイデル法SOR法の3手法の比較を行いました。

離散化式はそれぞれ似ていますが、3つの手法の中でSOR法が圧倒的にはやいということも確かめられ、さらに動画による違いも確認することができました。

ヤコビ法は並列計算ができ、その他はできないので一概にどの手法が優れているかは言えませんが、アルゴリズムで挙動が劇的に変わる様子を目の当たりにできたのは面白かったです。

(forループ使わない方法、筆者はまだヤコビ法しかかけない…。SOR法をforループ最低限版で書きたいですがまだ思いつきません。ダレカオシエテクダサイ…。)

参考文献

[1]:斎藤恭一, なっとくする偏微分方程式, 講談社, (2008), pp100 [2]:山崎郭滋, 偏微分方程式の数値解法入門, 森北出版,(2015), pp28

ラプラス方程式を数値計算できました!2次元は可視化が派手で面白いですね!
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

SNSでもご購読できます。

コメントを残す

*