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)

を用いています。

Page updated

Google Sites

Report abuse