Part.19
2次元拡散方程式（熱伝導方程式）の数値解法
空間２次元の熱伝導方程式を考えます。１次元の問題が、棒の熱分布を考えていたのに対し、２次元の問題では、薄い板の熱分布を考えていると想像してください。
問題は、
Find u(x,y,t) s.t.
となります。これを数値計算する為に離散化を行います。一次元の場合と同様なので、途中は省略すると、
という式が得られます。まとめると、
となります。特に、dx=dy として、
とかけば、
なる安定性条件を満たす必要があるので、注意が必要です（一次元の場合よりも、条件が厳しくなっています）。
これを解く為には、２次元配列が必要です。double 型の２次元配列は次のように宣言します。
double A[100][100];
これで、100x100の配列を作成できます。配列の要素は A[1][4] といったふうに参照できます。これで行列Aの１行４列目の要素を参照できます。（もちろん大きさには制限があります。あまりに巨大なものは扱えません。）
double A[101][101];
と宣言すると、A[0][0], A[0][1], ...., A[0][99],.....,A[100][100] と101x101の配列となります。
同様に、
double u[101][101];
double new_u[101][101];
という配列を作れば、一次元の熱方程式のアルゴリズムほぼそのままで計算できます。次の例では、D=1, Lx=Ly=1 の場合の例です。ただし完成版ではなく、太字のコメント部分が抜けています。
/* ２次元熱伝導方程式のアニメーション *//* D.Ueyama *//* 未完成版 */#include <glsc.h>#include <math.h>#include <stdio.h>/* 定数の設定 */#define PI (3.1415926)/* 空間刻数 */#define IMAX (100)#define JMAX (IMAX)/* TMAX */#define TMAX (0.05)/* Lx, Ly */#define Lx (1.0)#define Ly (Lx)/* 画面表示間隔 */#define INTV (100)/* 仮想座標系の下限・上限 */#define BOT (-1.0)#define TOP ( 1.0)int main(){ /* プログラム中で使用する変数の型宣言 */ int i, j, k; int KMAX; double u[IMAX+1][JMAX+1], new_u[IMAX+1][JMAX+1], t, x, y, dt, dx, dy, ax, ay; char text[256]; /* dt, dx 等を設定 */ dx = Lx/IMAX; dy = Ly/JMAX; /* ここでは dx = dy を仮定している */ /* ax が 1/6 となるようにしているが */ /* 実際には dx,dy のうち小さい方から dt を決めるべきである */ dt = (dx*dx)/6.0; ax = dt/(dx*dx); ay = dt/(dy*dy); KMAX = TMAX/dt; /* GLSCの初期設定 */ g_init("Graph", 200.0, 200.0); g_device(G_DISP); /* 仮想座標系を一つ用意 （今回は使わない）*/ g_def_scale(0, 0.0, Lx, BOT, TOP, 15.0, 10.0, 175.0, 80.0); g_cls(); g_sel_scale(0); /* 初期値設定 */ /* 境界条件設定 */ /* 時間方向へのループ（アニメーション）*/ for(k = 0; k <= KMAX; k++) { t = k*dt; /* INTV 間隔でグラフを表示 */ if( k%INTV == 0 ) { /* 一つ前の時間のグラフを消す */ g_cls(); /* テキストの表示 */ sprintf(text, "2D Heat Equation (%d*%d) by u11****", IMAX, JMAX); g_text(50.0,7.0, text); /* t の表示 */ sprintf(text, "t = %7.5lf", t); g_text(85.0,20.0, text); /* 配列 u をグラフとして描く （鳥瞰図）*/ g_line_color(G_RED); g_line_width(1.0); g_hidden(100.0, 100.0, 50.0, 0.0, 1.0, 500.0, 40.0, 35.0, 10.0, 10.0, 180.0, 180.0, (double *)u, IMAX+1, JMAX+1, 1, 0, 5, 5); if( k == 0 ) { /* 初期値を良くみたいので */ /* 初期値の場合は表示しクリックを待つ */ g_sleep(G_STOP); } /* 初期値以外の場合は 0.01 秒間止る */ /* これによりアニメーションのちらつきをある程度おさえる */ else { g_sleep(0.01); } } /* 空間方向へのループ */ /* 時間を1ステップ進める */ /* u から new_u を求める */ /* new_u を u にコピー */ } /* 計算終了後マウスがクリックされるまで停止 */ g_sleep(G_STOP);}課題１
初期値は境界条件を満たす適当なものを与えるものとし、上記プログラムを完成させよ。
課題２
プログラム中にコメントとして書いてあるが、ここでは dx = dy を仮定して ax = 1/6 となるように dt を決定している。このままだと、 dx > dy といった状況（例えば、 IMAX=100 で JMAX=200 とした場合など）では安定性条件を破る可能性がある（安定した計算のためには ax, ay それぞれが 1/4 より小さいことが必要）。dx と dy の大小関係から ax および ay のうち大きい方が 1/6 となるように dt を決定するプログラムに変更せよ。
上記プログラムでは図を鳥瞰図で表示しています。各引数の意味は、GLSCマニュアルを参照してください。適当に見る角度等は変えてもらって結構です。
g_hidden(100.0, 100.0, 50.0, 0.0, 1.0, 500.0, 40.0, 35.0, 10.0, 10.0, 180.0, 180.0,
(double *)u, IMAX+1, IMAX+1, 1, 0, 5, 5);
ただし、上記例は
g_init("Graph", 200.0, 200.0);
にて開いたウィンドウを仮定しています。
上の例にある g_hidden の引数中に (double *)u なるものがありますが、単に u と書いても正常に動きます。ただし、その場合コンパイル時に警告がでます。この (float *) なる記述はキャストと呼ばれるもので、変数の型をかえる場合に使います。今回の場合、 g_hidden 関数の第１３引数は float * 型を要求するのですが、ここでの u は２次元配列であるため、 float ** 型となり、キャストが必要となります。
２次元熱方程式のアニメーションサンプルを用意しました。このサンプルでは、初期値として、
sin(3*PI*x)*sin(2*PI*y)
を用いています。