おまけ回
2次元偏微分方程式の数値解法
（2022/03/15 最終修正）
空間2次元の熱伝導方程式を考える。1次元の問題が、棒の熱分布を考えていたのに対し、2次元の問題では、薄い板の熱分布を考えていると想像すればよい。
問題は、
Find u(x,y,t) s.t.
となる。これを数値計算する為に離散化を行う。1次元の場合と同様なので、途中は省略すると、
という式が得らる。まとめると、
となる。特に、dx=dy として、
とかけば、
なる安定性条件を満たす必要があるので、注意が必要（1次元の場合よりも、条件が厳しくなっている）。
これを解く為には、2次元配列が必要。さすがに NumPy を使おうと思う。NumPy を使うには、次の import 文が必要となる。
import numpy as np
その上で、
A = np.zeros((100, 100))
とすれば、100x100の配列を作成でき、変数の値は全て 0.0 となる。配列の要素は A[1][4] というふうに参照する。これで行列Aの1行4列目の要素を参照できる。（配列の添え字は0からはじまるので、0行3列目という言い方を許していただいた場合の言い方。それを許さないというのであれば、2行5列目の要素ということになる。）
A = np.zeros((101, 101))
と宣言すると、A[0][0], A[0][1], ...., A[0][99],.....,A[100][100] と101x101の配列となる。
さて、1次元問題に対して2次元問題は空間方向が一つ増えるので、繰り返し文が一つ深くなる。それだけといえばそれだけである。
次に、数値計算を Python で主に行っている人には笑われるだろうが、Gray-Scottモデルの空間2次元シミュレーションのコードを載せ、その解説をすることにする。Python で、こんなことをしたら遅くてたまらないですよ！と言われるやり方である。ただ、C言語などで書く場合との対応はつきやすいので、プログラミングの練習としては良いだろう。実際、最後の方でC言語で書き換えを行う。
Python 大好きで、もう Python で数値計算もします！という人は、NumPy, Numba というキーワードで検索すると色々情報が得られるので、各人試すとよい。それらを使えば高速化はできる（らしい）。
Code-gs2d.py : 2次元 Gray-Scott モデル（反応拡散系）のシミュレーション
# Code-gs2d.py
import numpy as np
import math
import matplotlib.pyplot as plt
import matplotlib.animation as animation
T = 1000.0
M = 2000
N = 100
Du = 1.0e-5
Dv = 2.0e-5
F = 0.02
kappa = 0.059
L = 1.0
Lc = L/2
h = L/float(N)
tau = T/float(M)
alphau = Du*tau/(h*h)
alphav = Dv*tau/(h*h)
INTV = 50
print(max(alphau, alphav))
def u0(x, y):
if x > Lc-0.1 and x < Lc+0.1 and y > Lc-0.1 and y < Lc+0.1 :
return 1/4.0
else:
return 0.0
def v0(x, y):
if x > Lc-0.1 and x < Lc+0.1 and y > Lc-0.1 and y < Lc+0.1 :
return 0.5
else:
return 1.0
def fu(U, V):
return (U*U*V - (F + kappa)*U)
def fv(U, V):
return (-U*U*V + F*(1.0 - V))
def setBC():
for i in range(1, N + 1):
u[i][0] = u[i][1]
u[i][N+1] = u[i][N]
v[i][0] = v[i][1]
v[i][N+1] = v[i][N]
for j in range(1, N + 1):
u[0][j] = u[1][j]
u[N+1][j] = u[N][j]
v[0][j] = v[1][j]
v[N+1][j] = v[N][j]
return
u = np.zeros((N+2, N+2))
new_u = np.zeros((N+2, N+2))
v = np.zeros((N+2, N+2))
new_v = np.zeros((N+2, N+2))
X, Y = np.meshgrid(np.linspace(0, L, N+2), np.linspace(0, L, N+2))
fig = plt.figure()
ims = []
for i in range(1, N + 1):
x = (i - 0.5)*h
for j in range(1, N + 1):
y = (j - 0.5)*h
u[i][j] = u0(x,y)
v[i][j] = v0(x,y)
setBC()
im = plt.contourf(X, Y, u, 100)
ims.append(im.collections)
for k in range(1, M+1):
print(k)
t = k*tau
setBC()
for i in range(1, N + 1):
for j in range(1, N + 1):
new_u[i][j] = u[i][j] + alphau*(u[i-1][j] + u[i+1][j] + u[i][j-1] + u[i][j+1] - 4*u[i][j]) + fu(u[i][j], v[i][j])*tau
new_v[i][j] = v[i][j] + alphav*(v[i-1][j] + v[i+1][j] + v[i][j-1] + v[i][j+1] - 4*v[i][j]) + fv(u[i][j], v[i][j])*tau
for i in range(1, N + 1):
for j in range(1, N + 1):
u[i][j] = new_u[i][j]
v[i][j] = new_v[i][j]
if (k + 1)%INTV == 0:
im = plt.contourf(X, Y, u, 100)
ims.append(im.collections)
plt.colorbar()
ani = animation.ArtistAnimation(fig, ims, interval=100, repeat=False)
#plt.show()
ani.save("gs2d.gif", writer="pillow", fps=10)
#ani.save("gs2d.mp4", writer="ffmpeg", fps=10)
動画中のコードには間違いがあり、それについては本ページの最後の方（2次元波動方程式についての前）に記述がある。上記コードは修正済み。
上の動画は、Code-gs2d.py 中の T, M をそれぞれ10倍してシミュレーションした結果（かなり時間を使うのと、画像をリスト内に保存するので、かなりのメモリを使用するので注意が必要）。長時間計算すると、2次元のこのパラメータの場合にはスポットの消滅生成が繰り返される。その現象については、以下の論文で検討した。
Spatio-temporal chaos for the Gray-Scott Model, Y.Nishiura and D.Ueyama, Physica D 150 (2001) , pp.137-162 .
C言語でのシミュレーションに移行してみる
工夫すればPythonでも速くなるようだが、やはり本格的な？数値シミュレーションをするならば、C言語への移行を考えるべきだと思う。問題は、数値シミュレーション結果の可視化であるが、そこをプログラムから分離し、数値データをファイルとして保存し、それをアプリを使って可視化するという方法をとることで、C言語への移行を実現してみようと思う。まずは、Pythonからmatlibplot のコードを削除し、ファイルに数値データを書きだす様に変更する。その後、そのコードをC言語に書き直す。PythonとC言語の関係性を理解していただければ、C言語への移行はそれほど大変では無い事がわかるだろう。データの可視化には ParaView を用いる。ダウンロードしておこう。macOS, Windows, Linux 用がある（ソースコードも公開されている）。
Python コードのスリム化（グラフィックス部分の消去）と VTK 形式での書きだしコードの追記
データファイルの書き出しは、適当に出力してから可視化アプリが読み込める形式に変更する方法と、書き出し時点で可視化アプリに合わせて書きだす方法があるだろう。ここでは、後者の立場を取り、ParaView が読み込むことができる VTK 形式でデータをファイルとして書きだす。各時間ステップ毎に連番でファイルを作成することで、ParaView でアニメーション化できる。実は VTK 形式に NumPy 配列を書きだすパッケージが幾つかあることを調べているうちに発見したが、他のパッケージにあまり頼らないことで、次のC言語への移行をスムーズに行いたいと思い、汚らしいコードだが VTK に書きだす関数を作成した。なんだ、こんなもんかと思うのも大事だろう。
Code-gs2d-vtk.py : matplotlib を使うのをやめて、VTK形式でデータを保存する。ParaViewを使ってみる。
# Code-gs2d-vtk.py
import numpy as np
import math
T = 10000.0
M = 20000
N = 100
Du = 1.0e-5
Dv = 2.0e-5
F = 0.02
kappa = 0.059
L = 1.0
Lc = L/2
h = L/float(N)
tau = T/float(M)
alphau = Du*tau/(h*h)
alphav = Dv*tau/(h*h)
INTV = 50
cc = 0
print(max(alphau, alphav))
def writeVtk(count):
fp = open("data_" + str(count) + ".vtk", "w")
fp.write("# vtk DataFile Version 4.1 \n")
fp.write("COMMENT\n")
fp.write("ASCII\n")
fp.write("DATASET STRUCTURED_POINTS \n")
fp.write("DIMENSIONS " + str(N) + " " + str(N) + " 1 \n")
fp.write("ORIGIN 0 0 0\n")
fp.write("SPACING " + str(h) + " " + str(h) + " 0 \n")
fp.write("POINT_DATA " + str(N*N) + "\n")
fp.write("SCALARS U double 1\n")
fp.write("LOOKUP_TABLE default\n")
for i in range(1, N + 1):
for j in range(1, N + 1):
fp.write(str(u[i][j]) + "\n")
fp.close()
def u0(x, y):
if x > Lc-0.1 and x < Lc+0.1 and y > Lc-0.1 and y < Lc+0.1 :
return 1/4.0
else:
return 0.0
def v0(x, y):
if x > Lc-0.1 and x < Lc+0.1 and y > Lc-0.1 and y < Lc+0.1 :
return 0.5
else:
return 1.0
def fu(U, V):
return (U*U*V - (F + kappa)*U)
def fv(U, V):
return (-U*U*V + F*(1.0 - V))
def setBC():
for i in range(1, N + 1):
u[i][0] = u[i][1]
u[i][N+1] = u[i][N]
v[i][0] = v[i][1]
v[i][N+1] = v[i][N]
for j in range(1, N + 1):
u[0][j] = u[1][j]
u[N+1][j] = u[N][j]
v[0][j] = v[1][j]
v[N+1][j] = v[N][j]
return
u = np.zeros((N+2, N+2))
new_u = np.zeros((N+2, N+2))
v = np.zeros((N+2, N+2))
new_v = np.zeros((N+2, N+2))
for i in range(1, N + 1):
x = (i - 0.5)*h
for j in range(1, N + 1):
y = (j - 0.5)*h
u[i][j] = u0(x,y)
v[i][j] = v0(x,y)
setBC()
for k in range(1, M+1):
t = k*tau
setBC()
for i in range(1, N + 1):
for j in range(1, N + 1):
new_u[i][j] = u[i][j] + alphau*(u[i-1][j] + u[i+1][j] + u[i][j-1] + u[i][j+1] - 4*u[i][j]) + fu(u[i][j], v[i][j])*tau
new_v[i][j] = v[i][j] + alphav*(v[i-1][j] + v[i+1][j] + v[i][j-1] + v[i][j+1] - 4*v[i][j]) + fv(u[i][j], v[i][j])*tau
for i in range(1, N + 1):
for j in range(1, N + 1):
u[i][j] = new_u[i][j]
v[i][j] = new_v[i][j]
if (k + 1)%INTV == 0:
print(cc)
writeVtk(cc)
cc = cc + 1
動画中のコードには間違いがあり、それについては本ページの最後の方（2次元波動方程式についての前）に記述がある。上記コードは修正済み。
Code-gs2d-vtk.c : PythonのコードをC言語で書き直す。Pythonで2000秒かかっていた計算が4秒になった。
//# Code-gs2d-vtk.py
#include <stdio.h>
#define T (10000.0)
#define M (20000)
#define N (100)
#define Du (1.0e-5)
#define Dv (2.0e-5)
#define F (0.02)
#define kappa (0.059)
#define L (1.0)
#define Lc (L/2)
#define h (L/(double)N)
#define tau (T/(double)M)
#define alphau (Du*tau/(h*h))
#define alphav (Dv*tau/(h*h))
#define INTV (50)
int cc = 0;
double u[N+2][N+2];
double new_u[N+2][N+2];
double v[N+2][N+2];
double new_v[N+2][N+2];
void writeVtk(int count) {
FILE *fp;
char fname[256];
sprintf(fname, "data_%d.vtk", count);
fp = fopen(fname, "w");
fprintf(fp, "# vtk DataFile Version 4.1 \n");
fprintf(fp, "COMMENT\n");
fprintf(fp, "ASCII\n");
fprintf(fp, "DATASET STRUCTURED_POINTS \n");
fprintf(fp, "DIMENSIONS %d %d 1\n", N, N);
fprintf(fp, "ORIGIN 0 0 0\n");
fprintf(fp, "SPACING %f %f 0 \n", h, h);
fprintf(fp, "POINT_DATA %d\n", N*N);
fprintf(fp, "SCALARS U double 1\n");
fprintf(fp, "LOOKUP_TABLE default\n");
for(int i = 1; i < N + 1; i++) {
for(int j = 1; j < N + 1; j++) {
fprintf(fp, "%f\n", u[i][j]);
}
}
fclose(fp);
}
double u0(double x, double y) {
if (x > Lc-0.1 && x < Lc+0.1 && y > Lc-0.1 && y < Lc+0.1) {
return 1/4.0;
}
else {
return 0.0;
}
}
double v0(double x, double y) {
if (x > Lc-0.1 && x < Lc+0.1 && y > Lc-0.1 && y < Lc+0.1) {
return 0.5;
}
else {
return 1.0;
}
}
double fu(double U, double V){
return (U*U*V - (F + kappa)*U);
}
double fv(double U, double V){
return (-U*U*V + F*(1.0 - V));
}
void setBC(){
for(int i = 1; i < N + 1; i++) {
u[i][0] = u[i][1];
u[i][N+1] = u[i][N];
v[i][0] = v[i][1];
v[i][N+1] = v[i][N];
}
for(int j = 1; j < N + 1; j++) {
u[0][j] = u[1][j];
u[N+1][j] = u[N][j];
v[0][j] = v[1][j];
v[N+1][j] = v[N][j];
}
}
int main() {
double x, y, t;
for(int i = 1; i < N + 1; i++) {
x = (i - 0.5)*h;
for(int j = 1; j < N + 1; j++) {
y = (j - 0.5)*h;
u[i][j] = u0(x,y);
v[i][j] = v0(x,y);
}
}
setBC();
for(int k = 1; k < M + 1; k++) {
t = k*tau;
setBC();
for(int i = 1; i < N + 1; i++) {
for(int j = 1; j < N + 1; j++) {
new_u[i][j] = u[i][j] + alphau*(u[i-1][j] + u[i+1][j] + u[i][j-1] + u[i][j+1] - 4*u[i][j]) + fu(u[i][j], v[i][j])*tau;
new_v[i][j] = v[i][j] + alphav*(v[i-1][j] + v[i+1][j] + v[i][j-1] + v[i][j+1] - 4*v[i][j]) + fv(u[i][j], v[i][j])*tau;
}
}
for(int i = 1; i < N + 1; i++) {
for(int j = 1; j < N + 1; j++) {
u[i][j] = new_u[i][j];
v[i][j] = new_v[i][j];
}
}
if((k + 1)%INTV == 0) {
printf("%d\n", cc);
writeVtk(cc);
cc = cc + 1;
}
}
}
C言語への書き換えを終えて、結果がPythonのものと異なる事に気が付きましたが、Pythonプログラムの字下げが一カ所間違っていたことによるものでした。動画中では倍精度がなんたらと言っていますが、そうではありません。動画はそのままにしておきますが、本ページ内の2つのPythonプログラムは書き換えました（太字の return 文）。字下げが一つ深かったために、for 文を完結する事無く途中で関数を抜けていました。よって、境界条件が正しく設定されておらず対称性のズレが生じていました。それでもおおよそ正しい結果になる所が数値計算の怖いところではあります。ご注意下さい。
openprocessing.org に、同様のシミュレーションのコードを置いている。次にそれらを貼っておく。こちらは、JavaScript で書かれている。マウスでクリックすると、その周りのスポットを消す事ができる。</>を押すとコードをみる事ができる。パラメータの変更もできたと思う。
2次元波動方程式について
こちら Python コードを載せておく。C言語への書き換えは、先の反応拡散系の例を見てご自身でやってみてほしい。良い練習になるだろう。VTKファイルを ParaView で可視化したものを載せておく。
Code-wave2d-vtk.py : 2次元波動方程式
# Code-wave2d-vtk.py
import numpy as np
import math
T = 3.0
M = 3000
tau = T/M
L = 2.0
N = 200
h = L/N
C = 1.0
Lambda = C*tau/h
INTV = 10
cc = 0
print(Lambda)
def writeVtk(count):
fp = open("data_" + str(count) + ".vtk", "w")
fp.write("# vtk DataFile Version 4.1 \n")
fp.write("COMMENT\n")
fp.write("ASCII\n")
fp.write("DATASET STRUCTURED_POINTS \n")
fp.write("DIMENSIONS " + str(N) + " " + str(N) + " 1 \n")
fp.write("ORIGIN 0 0 0\n")
fp.write("SPACING " + str(h) + " " + str(h) + " 0 \n")
fp.write("POINT_DATA " + str(N*N) + "\n")
fp.write("SCALARS U double 1\n")
fp.write("LOOKUP_TABLE default\n")
for i in range(1, N + 1):
for j in range(1, N + 1):
fp.write(str(old_u[i][j]) + "\n")
fp.close()
def f(x, y):
return math.exp(-100*((x-0.8)*(x-0.8) + (y-0.5)*(y-0.5)))
def g(x, y):
return 0
x, y = np.array([0.0]*(N+2)), np.array([0.0]*(N+2))
u, old_u, new_u = np.zeros((N+2, N+2)), np.zeros((N+2, N+2)), np.zeros((N+2, N+2))
#初期値設定
for j in range(1, N+1):
y[j] = (j - 0.5)*h
for i in range(1, N+1):
x[i] = (i - 0.5)*h
old_u[i][j] = f(x[i], y[j])
#境界条件
for j in range(1, N+1):
old_u[0][j] = old_u[1][j]
old_u[N+1][j] = old_u[N][j]
for i in range(1, N+1):
old_u[i][0] = old_u[i][1]
old_u[i][N+1] = old_u[i][N]
for j in range(1, N+1):
for i in range(1, N+1):
u[i][j] = old_u[i][j] + tau*g(x[i], y[j]) + (Lambda*Lambda/2)*(old_u[i-1][j] - 2*old_u[i][j] + old_u[i+1][j]) + (Lambda*Lambda/2)*(old_u[i][j-1] - 2*old_u[i][j] + old_u[i][j+1])
t = 0
for k in range(0, M+1):
t = k*tau
#境界条件
for j in range(1, N+1):
u[0][j] = u[1][j]
u[N+1][j] = u[N][j]
for i in range(1, N+1):
u[i][0] = u[i][1]
u[i][N+1] = u[i][N]
for j in range(1, N+1):
for i in range(1, N+1):
new_u[i][j] = 2*u[i][j] - old_u[i][j] + (Lambda*Lambda)*(u[i-1][j] - 2*u[i][j] + u[i+1][j]) + (Lambda*Lambda)*(u[i][j-1] - 2*u[i][j] + u[i][j+1])
for j in range(1, N+1):
for i in range(1, N+1):
old_u[i][j] = u[i][j]
u[i][j] = new_u[i][j]
if k%INTV == 0:
print(cc)
writeVtk(cc)
cc = cc + 1