こんにちは(@t_kun_kamakiri)。
前回の記事では、「2次元のバーガース方程式」をGoogle Colabでアニメーション作成を行いました。
前回の記事はこちら
今回は、以下のような２次元のラプラス方程式をPythonで実装します。
なぜ、ラプラス方程式の解法を知る必要があるのかというと、ナビエストークス方程式を解く際に必要になるからです。（正確には、必要なのはポアソン方程式です。）
本記事では、そちらも詳しく解説します。
本件の基本的な内容はこちらのサイトにそってやっていきます。
この記事ではこんな人を対象にしています。
こんな人が対象
- Pythonを使い始めたけどどう使うかわからない
- 流体の数値計算をはじめて勉強する人
本記事シリーズの最終目標は、ナビエストークス方程式をPythonで実装することですが、いきなりナビエストークス方程式を実装するのは難しいので各項の意味を確認しながら進めていきたいと思います。
今回の内容はこちら
- ２次元のラプラス方程式の数値的解法をなぜ知る必要があるのか
- ２次元のラプラス方程式の数値計算をPythonで実装
カマキリ
今日もゆる～く学びます(^^)/
2次元のラプラス方程式
今回は以下のラプラス方程式ついてPythonを用いて実装して数値計算を行いたいと思います。
\begin{align*}
\frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = 0
\tag{1}
\end{align*}
\frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = 0
\tag{1}
\end{align*}
いきなりラプラス方程式の話をされても、なぜここでラプラス方程式の話をするのかと、疑問に思う人もいることでしょう！
- ラプラス方程式とポアソン方程式の違い
- ナビエストークス方程式の数値計算にポアソン方程式を解く場面がある
ラプラス方程式とポアソン方程式の違い
ラプラス方程式の形をしていても、一般的には右辺が0ではない場合、
\begin{align*}
\frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = f(x,y)\tag{2}
\end{align*}
\frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = f(x,y)\tag{2}
\end{align*}
をポアソン方程式と呼ばれています。
ポアソン方程式として最も有名な式としては、「電荷の分布\(\rho(\boldsymbol{r})\)があります。
周りの静電場\(\boldsymbol{E}\)」マクスウェル方程式\(\nabla\cdot\boldsymbol{E}=\frac{\rho(\boldsymbol{r})}{\epsilon}\)を満たし、スカラー場\(\phi\)としたとき、スカラー場の勾配が電場\(\boldsymbol{E}=-\nabla\phi\)なので、
\begin{align*}
\nabla^2 \phi =-\frac{\rho(\boldsymbol{r})}{\epsilon}
\end{align*}
\nabla^2 \phi =-\frac{\rho(\boldsymbol{r})}{\epsilon}
\end{align*}
電荷が0の場合は、ラプラス方程式になるというわけですね。
明らかにポアソン方程式の難しそうなので、まずは理解のために簡単そうなラプラス方程式について考えてみましょう。
というのが本記事の狙いでもあります。
ナビエストークス方程式の数値計算にポアソン方程式を解く場面がある
ところで、なぜラプラス方程式をここで持ち出したかというと、
ナビエストークス方程式の数値計算にポアソン方程式を解く場面があるからです。
非圧縮性流体として現象を取り扱う場合には、基本的には以下の2つの式を連立して解きます。
実際の流れは、すべての現象が(3)(4)式というわけではないで。
「ニュートン流体」「非圧縮の流れ（近似）」としているとても限定的な流れであることを理解しておきましょう。
もっと詳しく勉強したい方は☟こちらの記事を参考にしてください。
もっと詳しく勉強したい方は☟こちらの記事を参考にしてください。
前回の記事はこちら
今回は、2次元で計算を行うため少々くどいですが、成分ごとの式を書き下すことにします。
- ナビエストークスからは2本の式
- 連続の式から1つの式
が出てくるので、合計で3本の式があります。
未知数は、\(\boldsymbol{v}\)の2つと圧力\(p\)の合計3つです。
なので、原理的に未知数である物理量は全て求まることになります。
それを各時刻で求めていけば良いという事になります。
しかし、言葉で簡単に言っても、特殊な境界条件や近似を用いないとナビエストークス方程式は解析的に取り扱うことができないのです。
つまり、ナビエストークス方程式には一般解が見つかっていないので、実際には数値計算を使って偏微分方程式を取り扱うしかないのが現状です。
しかし、連立方程式を解けばいいと言われてもちょっと困った・・・
なぜなら、(5)(6)式で求めた流速は(7)式を満たす必要があります。
そのように流速を選びながら圧力を同時に決めなければなりません。
考えただけでもちょっと工夫が必要ですよね。。
そこで、以下のようにして考えを改めます。
圧力に関しては非圧縮条件を満たすようにポアソン方程式を解きつつ、求めた圧力に関してナビエストークス方程式を解き流速を求めるという手法をとります。
言葉で説明してもよくわからないかもしれないので、実際にやります(^^)/
(5)式を\(y\)で偏微分
\begin{align*}
&\frac{\partial }{\partial x}\frac{\partial u}{\partial t}+\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+u\frac{\partial^2 u}{\partial x^2}+ \frac{\partial v}{\partial x}\frac{\partial u}{\partial y}+v\frac{\partial^2 u}{\partial x\partial y} \\
&= -\frac{1}{\rho}\frac{\partial^2 p}{\partial x^2}+\nu \frac{\partial }{\partial x}\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2} \right)\tag{8}
\end{align*}
&\frac{\partial }{\partial x}\frac{\partial u}{\partial t}+\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+u\frac{\partial^2 u}{\partial x^2}+ \frac{\partial v}{\partial x}\frac{\partial u}{\partial y}+v\frac{\partial^2 u}{\partial x\partial y} \\
&= -\frac{1}{\rho}\frac{\partial^2 p}{\partial x^2}+\nu \frac{\partial }{\partial x}\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2} \right)\tag{8}
\end{align*}
(6)式を\(x\)で偏微分
\begin{align*}
&\frac{\partial }{\partial y}\frac{\partial v}{\partial t}+\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+u\frac{\partial^2 v}{\partial y\partial x}+ \frac{\partial v}{\partial y}\frac{\partial v}{\partial y}+v\frac{\partial^2 v}{\partial x\partial y^2} \\
&=-\frac{1}{\rho}\frac{\partial^2 p}{\partial y^2}+\nu \frac{\partial }{\partial y}\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2} \right)\tag{9}
\end{align*}
&\frac{\partial }{\partial y}\frac{\partial v}{\partial t}+\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+u\frac{\partial^2 v}{\partial y\partial x}+ \frac{\partial v}{\partial y}\frac{\partial v}{\partial y}+v\frac{\partial^2 v}{\partial x\partial y^2} \\
&=-\frac{1}{\rho}\frac{\partial^2 p}{\partial y^2}+\nu \frac{\partial }{\partial y}\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2} \right)\tag{9}
\end{align*}
このようにして辺々足して、非圧縮条件を使えば、
\begin{align*}
\frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2} = -\rho\left(\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+2\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+\frac{\partial v}{\partial y}\frac{\partial v}{\partial y} \right)\tag{10}
\end{align*}
\frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2} = -\rho\left(\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+2\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+\frac{\partial v}{\partial y}\frac{\partial v}{\partial y} \right)\tag{10}
\end{align*}
計算はだるいですが、一度やってみてください。
つまり、(5)(6)(7)を解く代わりに以下の式を解けばよいという事になります。
【2次元での計算】※実際に解く方程式
ナビエストークス方程式
\begin{align*}
\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial x}+\nu \left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2} \right) \tag{5}
\end{align*}
\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial x}+\nu \left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2} \right) \tag{5}
\end{align*}
\begin{align*}
\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial y}+\nu\left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}\right)\tag{6}
\end{align*}
\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial y}+\nu\left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}\right)\tag{6}
\end{align*}
ポアソン方程式※連続の式（非圧縮条件）
\begin{align*}
\frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2} = -\rho\left(\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+2\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+\frac{\partial v}{\partial y}\frac{\partial v}{\partial y} \right)\tag{11}
\end{align*}
\frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2} = -\rho\left(\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+2\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+\frac{\partial v}{\partial y}\frac{\partial v}{\partial y} \right)\tag{11}
\end{align*}
この3つの式を解けばよいという事になります。
で、3つ目の式がポアソン方程式です。
ポアソン方程式は数値的な解法としては有名なので、このような形に式変形しました。
速度にいては、ナビエストークス方程式を解きつつ、圧力に対してポアソン方程式を連立させながら解くという事になります。
だから、ここでポアソン方程式を学んでおこうというわけです。
ですが、
ポアソン方程式よりも簡単な右辺が0であるラプラス方程式を学びましょうというのが本記事の趣旨です。
ラプラス方程式の離散化
くどいですが、今回の記事では☟こちらの方程式を扱う・・・・
ポアソン方程式※連続の式（非圧縮条件）
\begin{align*}
\frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2} = b(x,y)\tag{12}
\end{align*}
\frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2} = b(x,y)\tag{12}
\end{align*}
\begin{align*}
b(x,y)=-\rho\left(\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+2\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+\frac{\partial v}{\partial y}\frac{\partial v}{\partial y} \right)\tag{13}
\end{align*}
b(x,y)=-\rho\left(\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+2\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+\frac{\partial v}{\partial y}\frac{\partial v}{\partial y} \right)\tag{13}
\end{align*}
わけではなく、\(b(x,y)=0\)とした場合に、
ラプラス方程式
\begin{align*}
\frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2} =0\tag{14}
\end{align*}
\end{align*}
これを空間微分に対して、2次精度で離散化すると、
\begin{align*}
\frac{p_{i+1, j}^n – 2p_{i,j}^n + p_{i-1,j}^n}{\Delta x^2} + \frac{p_{i,j+1}^n – 2p_{i,j}^n + p_{i, j-1}^n}{\Delta y^2} = 0
\end{align*}
\frac{p_{i+1, j}^n – 2p_{i,j}^n + p_{i-1,j}^n}{\Delta x^2} + \frac{p_{i,j+1}^n – 2p_{i,j}^n + p_{i, j-1}^n}{\Delta y^2} = 0
\end{align*}
となります。
つまり、☟これを解けばよいという事です。
\begin{align*}
p_{i,j}^n = \frac{\Delta y^2(p_{i+1,j}^n+p_{i-1,j}^n)+\Delta x^2(p_{i,j+1}^n + p_{i,j-1}^n)}{2(\Delta x^2 + \Delta y^2)}\tag{15}
\end{align*}
p_{i,j}^n = \frac{\Delta y^2(p_{i+1,j}^n+p_{i-1,j}^n)+\Delta x^2(p_{i,j+1}^n + p_{i,j-1}^n)}{2(\Delta x^2 + \Delta y^2)}\tag{15}
\end{align*}
もし、\(\Delta x = \Delta y\)とした場合は、
\begin{align*}
p_{i,j}^n = \frac{1}{4}\big(p_{i+1,j}^n+p_{i-1,j}^n+p_{i,j+1}^n + p_{i,j-1}^n\big)\tag{16}
\end{align*}
p_{i,j}^n = \frac{1}{4}\big(p_{i+1,j}^n+p_{i-1,j}^n+p_{i,j+1}^n + p_{i,j-1}^n\big)\tag{16}
\end{align*}
(16)式からわかるように、\((i,j)\)の格子点での圧力\(p_{i,j}\)はまわり4点の圧力の平均であることがわかります。
式を見ればわかる通り、ラプラス方程式は時間微分の項がありません。離散化した場合では、\(p^{n+1}\)がないってことです。
つまり、ラプラス方程式は平衡状態に達した場合における偏微分方程式です。
有限の空間で行われる数値計算においては、ラプラス方程式は境界条件によって空間分布の状態が決定するため「境界値問題」とも言われています。
要するに境界条件を決めたら物理量の空間分布は決まりますよってことです。
境界条件は、基本的には計算領域の一番端をどのようにするかで決まります。
今回行う境界条件は以下のように設定します。
【ディリクレ境界条件】
- \(p=0\) at
- \(x=0\)
- \(p=y\) at \(x=2\)
【ノイマン境界条件】
\(\frac{\partial p}{\partial y}=0\) at \(y=0, 1\)
境界条件を可視化すると☟このようになります。
この境界条件をもとにラプラス方程式を解くと、
\begin{align*}
p(x,y)=\frac{x}{4}-4\sum_{n=1,odd}^{\infty}\frac{1}{(n\pi)^2\sinh2n\pi}\sinh n\pi x\cos n\pi y\tag{17}
\end{align*}
p(x,y)=\frac{x}{4}-4\sum_{n=1,odd}^{\infty}\frac{1}{(n\pi)^2\sinh2n\pi}\sinh n\pi x\cos n\pi y\tag{17}
\end{align*}
となるとこちらのサイトでは言っていますが、
カマキリ
どうやったら(17)式になるんだろう・・・・って考え中です。
ラプラス方程式のこいつがある意味がわかりません_(._.)_ pic.twitter.com/jS0F7vzF8Y
— カマキリ🐲@物理ブログ書いている (@t_kun_kamakiri) June 17, 2020
Twitterでも質問してみましたが、情報を得ることができず・・・・誰かわかったら教えてください。
（ラプラス方程式は変数分離をして、境界条件を課すことでできるというのが僕の中の定石ですが、上手くいきませんでした_(._.)_）
（ラプラス方程式は変数分離をして、境界条件を課すことでできるというのが僕の中の定石ですが、上手くいきませんでした_(._.)_）
また別途、ラプラス方程式の解析解についてと数値計算の解法についての解説記事を買いたいと思います。
では、以下で実際にラプラス方程式をGoogle Colaboratoryを使ってPythonのコードを書きながら理解を深めていきたいと思います。
ラプラス方程式をPythonで実装
ひとつひとつ解説をしながらコードを書いていきます。
※関数にまとめておいた方がいい部分はあとでまとめたいと思います。
【ポイント】
- 必要なライブラリをインポートする。
コードを以下のように書きます。
1 2 3 4 | import numpy as np from matplotlib import pyplot as plt from matplotlib import pyplot, cm from mpl_toolkits.mplot3d import Axes3D |
NumpyとMatplotlibをインポートします。
記述ミスがなければエラーなく進みます。
では、次にラプラス方程式における「初期の状態」と「境界条件」を課します。
【ポイント】
- 初期状態を設定
- 境界条件を設定
コードを以下のように書きます。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 | nx = 31 ny = 31 c = 1 dx = 2 / (nx - 1) dy = 2 / (ny - 1) ##initial conditions p = np.zeros((ny, nx)) # create a XxY vector of 0's ##plotting aids x = np.linspace(0, 2, nx) y = np.linspace(0, 1, ny) #境界条件 p[:, 0] = 0 # p = 0 @ x = 0 p[:, -1] = y # p = y @ x = 2 p[0, :] = p[1, :] # dp/dy = 0 @ y = 0 p[-1, :] = p[-2, :] # dp/dy = 0 @ y = 1 |
これで、初期状態の圧力のプロファイルができています。
初期状態と言っても、境界条件を課しただけですので初期の状態には何の意味もありません。
あくまで境界条件をもとにラプラス方程式を解いて得られた圧力分布に意味があるので、初期状態はその最終状態を得るための種を与えているにすぎません。
そうは言っても数値計算では、初期の形状は確認しておいた方が良いですよね。
【ポイント】
Matplotlibを使って結果を可視化する。
コードを以下のように書きます。
1 2 3 4 5 6 7 8 9 10 11 | fig = plt.figure(figsize=(11, 7), dpi=100) ax = fig.gca(projection='3d') X, Y = np.meshgrid(x, y) surf = ax.plot_surface(X, Y, p[:], rstride=1, cstride=1, cmap=cm.viridis,linewidth=0, antialiased=False) ax.set_xlim(0, 2) ax.set_ylim(0, 1) ax.view_init(30, 225) ax.set_xlabel('$x$') ax.set_ylabel('$y$') ax.set_zlabel('$p$') |
【結果】
👆自分が想定していた境界条件になっていればOKです。
では、いよいよ上記で設定した境界条件をもとにラプラス方程式を解いてみましょう！
【ポイント】
\(p_{i,j}^n = \frac{1}{4}\big(p_{i+1,j}^n+p_{i-1,j}^n+p_{i,j+1}^n + p_{i,j-1}^n\big)\)について解くコードを実装します。
コードを以下のように書きます。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 | l1norm_target = 1e-4 l1norm = 1 pn = np.empty_like(p) while l1norm > l1norm_target: pn = p.copy() p[1:-1, 1:-1] = ((dy**2 * (pn[1:-1, 2:] + pn[1:-1, 0:-2]) + dx**2 * (pn[2:, 1:-1] + pn[0:-2, 1:-1])) / (2 * (dx**2 + dy**2))) p[:, 0] = 0 # p = 0 @ x = 0 p[:, -1] = y # p = y @ x = 2 p[0, :] = p[1, :] # dp/dy = 0 @ y = 0 p[-1, :] = p[-2, :] # dp/dy = 0 @ y = 1 l1norm = (np.sum(np.abs(p[:]) - np.abs(pn[:])) /np.sum(np.abs(pn[:]))) |
これで最終状態の圧力のプロファイルが完成しています。
何を行っているのかを「4×４」の格子状の点で簡単にして、回折します。
\begin{align*}
p_{i,j}^n = \frac{1}{4}\big(p_{i+1,j}^n+p_{i-1,j}^n+p_{i,j+1}^n + p_{i,j-1}^n\big)\tag{15}
\end{align*}
p_{i,j}^n = \frac{1}{4}\big(p_{i+1,j}^n+p_{i-1,j}^n+p_{i,j+1}^n + p_{i,j-1}^n\big)\tag{15}
\end{align*}
ですから、例えば(1,1)の格子点上の圧力は上記のように以下のようになります。
青文字は境界条件で値が与えられているので、既知の値です。
これを見ればわかるように、(i,j)の格子点での圧力はその周りの圧力の平均値として計算しています。
これを全ての格子点上で行うことを考えます。
今求めたいのは黒文字で書かれた圧力ですので、行列計算で求めることを考えます。
これで、未知の量について解くことができるのですが、
格子点数が多いとこのように行列計算で求めることは結構厳しいですよね。
この行列計算を行ってくれる方法としていろいろあるのですが、有名なのをいくつか挙げておきます。
- ガウス消去法（直接法）
- ガウス・ジョルダン法（直接法）
- ヤコビ法（反復法）
- ガウス・ザイデル法（反復法）
- SOR法（反復法）
などがあります。
詳しく勉強したい方は、以下の参考書がとてもわかりやすいです。
数値計算の闇について詳しく知りたい方は、☟こちらの参考書も合わせてお読みください。
その中でも今回は、反復法の中でも簡単な「ヤコビ法」で解くことにします。
収束判定は、
1 | l1norm = (np.sum(np.abs(p[:]) - np.abs(pn[:])) /np.sum(np.abs(pn[:]))) |
としています。
全ての格子点上の値の合計と前のステップでの値の合計とを比較して、「1e-4（0.01%)」以下の誤差になれば収束したと判定しています。
では、次に収束した圧力分布を可視化してみてみましょう。
【ポイント】
Matplotlibを使って結果を可視化する。
コードを以下のように書きます。
1 2 3 4 5 6 7 8 9 10 11 | fig = plt.figure(figsize=(11, 7), dpi=100) ax = fig.gca(projection='3d') X, Y = np.meshgrid(x, y) surf = ax.plot_surface(X, Y, p[:], rstride=1, cstride=1, cmap=cm.viridis,linewidth=0, antialiased=False) ax.set_xlim(0, 2) ax.set_ylim(0, 1) ax.view_init(30, 225) ax.set_xlabel('$x$') ax.set_ylabel('$y$') ax.set_zlabel('$p$') |
【結果】
このような結果になりました(^^)/
関数にまとめる
今のままで良いのですが、関数にまとめてコードをきれいにまとめておきます。
☟初期状態を設定します。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 | import numpy as np from matplotlib import pyplot as plt from matplotlib import pyplot, cm from mpl_toolkits.mplot3d import Axes3D nx = 31 ny = 31 c = 1 dx = 2 / (nx - 1) dy = 2 / (ny - 1) ##initial conditions p = np.zeros((ny, nx)) # create a XxY vector of 0's ##plotting aids x = np.linspace(0, 2, nx) y = np.linspace(0, 1, ny) #境界条件 p[:, 0] = 0 # p = 0 @ x = 0 p[:, -1] = y # p = y @ x = 2 p[0, :] = p[1, :] # dp/dy = 0 @ y = 0 p[-1, :] = p[-2, :] # dp/dy = 0 @ y = 1 |
次にラプラス方程式を解く部分を関数にまとめます。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 | def laplace2d(p, y, dx, dy, l1norm_target): l1norm = 1 pn = np.empty_like(p) i=0 while l1norm > l1norm_target: pn = p.copy() p[1:-1, 1:-1] = ((dy**2 * (pn[1:-1, 2:] + pn[1:-1, 0:-2]) + dx**2 * (pn[2:, 1:-1] + pn[0:-2, 1:-1])) / (2 * (dx**2 + dy**2))) p[:, 0] = 0 # p = 0 @ x = 0 p[:, -1] = y # p = y @ x = 2 p[0, :] = p[1, :] # dp/dy = 0 @ y = 0 p[-1, :] = p[-2, :] # dp/dy = 0 @ y = 1 l1norm = (np.sum(np.abs(p[:]) - np.abs(pn[:])) /np.sum(np.abs(pn[:]))) i +=1 print(i) return p |
ラプラス方程式を解いた後に、「return p」でpを返す関数を作成しました。
次に、結果を可視化する関数を用意します。
1 2 3 4 5 6 7 8 9 10 11 12 | def plot2D(x, y, p): fig = plt.figure(figsize=(11, 7), dpi=100) ax = fig.gca(projection='3d') X, Y = np.meshgrid(x, y) surf = ax.plot_surface(X, Y, p[:], rstride=1, cstride=1, cmap=cm.viridis,linewidth=0, antialiased=False) ax.set_xlim(0, 2) ax.set_ylim(0, 1) ax.view_init(30, 225) ax.set_xlabel('$x$') ax.set_ylabel('$y$') ax.set_zlabel('$p$') |
これで準備完成です(^^)/
初期の分布を見てみましょう。
1 | plot2D(x, y, p) |
【結果】
これは境界条件だけを可視化したもので、特に意味があるものではありませんね。
このあと、ラプラス方程式を解くことによって意味のある分布になります。
1 | p = laplace2d(p, y, dx, dy, 1e-4) |
これでラプラス方程式を解いたことになります。
いちおう収束のために何回ループしたのかも計算しておきました。
「1132回」収束のためにループ計算をしています。
では、最終状態の結果を見てみます。
1 | plot2D(x, y, p) |
【結果】
まとめ
今回は、2次元のラプラス方程式をPythonで実装してみました。
本記事のシリーズを読むことでPythonを使ったナビエストークス方程式の実装まではいけると思いますので、是非最後までお付き合いくださいm(__)m
次回の記事はこちら
おすすめの書籍紹介
ネット上には無料で多くのことが学べるのですが、ネット情報だけだと情報が偏ったりします。
やっぱり書籍を手に取って体系だって学んだ方が良い場合も多いです。
初心者でもので体系だって学ぶことができる、且つあまり難しすぎない書籍をここで紹介しておきます。
Pythonの完全初心者は書籍で学ぶとよい
Python自体が不安だという方は、こちらの書籍から勉強しても良いかも知れません。
☟こちらの本がPython初心者が挫折することなく勉強できる本です。
（本記事のようのPython使用環境と異なりますが、とてもわかりやすいので全く問題ありません）
僕自身もこの書籍ではじめは勉強しました。
小難しい余計なことが書いていないし、年末の一週間くらいで読むことができたので、初心者の方にはとてもお勧めです。
流体の勉強をしたい方
流体の数値計算以前に、流体力学を勉強しないと！！って思っている方は以下の参考書からはじめてみるのが良いと思います。
工学よりに書かれていて実用的な内容も触れながら、流体の理論的な内容も易しく書かれていますので大学１年生で偏微分を学んだ方にとってはちょうどよい参考書です。
数値流体の書籍
本記事で紹介している数値計算手法は、多くの数値解法のごく一部です。
特に、非線形方程式で記述できる流体現象のようなものは注意深く数値解法を選択しなくてはいけません。
以下の書籍は、値段が高いですが有限体積法をメインに数値解法がとても詳しく解説されています。
Follow @t_kun_kamakiri