流線関数を求める

前節で、導出した方程式を離散化してコンピューターに計算させる。まず、方程式のうち、流線関数と渦度の関係式(8 8)をMatlabを用いて解くことを考えよう。この式は、右辺の渦度、$\zeta$を既知として、左辺の流線関数、$\psi$を求めるものである。左辺は、流線関数のx、yについての2階の偏微分の和であり、この様なかたちの式を楕� 型方程式と呼ぶ。楕� 型の式では、境界条件が与えられれば、一つの解が得られる。式を差分化して次の様な離散方程式を得ることができる。離散の仕方には、様々な方法を考えることができるが、ここでは、単純に3点の計算点を用いて式を離散化してみる。

$$\frac{\psi (i-1,j)-2\psi(i,j)+\psi(i+1,j)}{\delta x^2} + \frac{\psi (i,j-1)-2\psi(i,j)+\psi(i,j+1)}{\delta z^2}=$$

$$\frac{u (i,j+1)-u(i,j-1)}{2\delta z} - \frac{w(i+1,j)-w(i-1,j)}{2\delta x},$$ (18)

前述の通り右辺の渦度は、既知として与える。この様な、楕� 型の方程式を解く為に通常のコンピューター言語では、 まず、流線関数の推定値を準備し、左辺を計算して右辺を 求める。求めた右辺と実際の渦度を比較して、その差が非常に小さくなるまで繰り返して、与えた渦度に対する流線関数を計算する。しかし ながら、Matlabは繰り返しの計算にはあまり向いていない。一方で、行列を直接計算できる。以下の様に式を行列で聞。

$$\left[ \begin{array}{ccccccccc} \frac{1}{\delta x^2+\delta y^2} &... ...)\\ :\\ :\\ :\\ \psi(N,M-2)\\ \psi(N,M-1)\\ \psi(N,M) \end{array}\right]=$$

$$\left[ \begin{array}{c} -\zeta(1,1)\\ -\zeta(1,2)\\ -\zeta(1,3)... ...:\\ :\\ :\\ -\zeta(N,M-2)\\ -\zeta(N,M-1)\\ -\zeta(N,M) \end{array}\right]$$ (19)

これは、行列の$Ax=b$のかたちをしており、N$\times$M個の未知数、$\psi$に対して、N$\times$M個の連立一次方程式を 雰、これを解く事で流線関数を求めることができる。前章で演習した、連立一次方程式を解く要領で解を求められそうだ。 しかしながら、左辺の行列は、N$\times$M行N$\times$M列の正方行列で、非常に大きい。しかも各行は、4から5つのゼロでない要素を持つが、それ以外は 全てゼロである。この様な行列をスパース行列と呼び、正方行列のほとんどの要素がゼロであるので、ゼロでない要素の情報だけが有用であり、 ゼロを含めて行列を準備するのは、非常に無駄である。Matlabでは、このようなスパース行列を効率よく格納し、それがつくる連立一次方程式を 解くためのツールが準備されている。そのうちの一つが、spconvertと呼ばれるMatlab関数である。これは、ゼロ以外の要素のみを、その行列における 位置情報とともにスパース行列として格納する機避持つ。このspconvertを用いて与えられた渦度から、流線関数を求めるMatlab プログラムは、以下の様なものになる。まず、正方行列を準備する。正方行列の要素は、格子間隔の関数となっており、 ここに流線関数、$\psi$の境界条件を組み込むことができる。例えば、x方向の$\psi$に関する離散式は、必ず$\psi(i+1,j)$や$\psi(i-1,j)$を含んでおり、 $i$を1-Nまで変化させた場合、格子外の点の情報が必要となる。ここに、周期的境界条件、 $\psi(0,j)=\psi(N,j)$や、 $\psi(N+1,j)=\psi(1,j)$を用いて、 存在しない点を存在する点で併ば、流線関数は上流と下流で同一の値を示すはずだ。同様に、z方向の境界条件$\psi=0$についても存在しない点の部分、 $\psi(i,0), \psi(i,M+1)$に関して、考慮しなければ、$\psi=0$を案に指定することができよう(該当する要素を0にする)。以下のMatlab Functionでは、この境界条件を考慮しつつ、 正方行列のゼロでない要素をA-Eとして生成し、それらが、行列の何行何列目に配置されるかを、例えばAに関しては、iA行、jA列として定義している。最終的に、何行何列という情報とゼロでない要素の情報をDATAというN$\times$M行、3列の行列として定義し、spconvertに渡してスパース行列を 効率よく格納している。

次に、上記の関数で準備した正方行列、MATRIXを用いて、式の右辺に相当する渦度を与えて、 $Ax=b$を解く。今、渦度、$\zeta$ (zeta)が$\psi$ (psi)と同様な格子状に与えられている場合、流線関数、$\psi$は 以下の様な、Matlab関数で得ることができよう。