3.5.2 熱方程式の初期値境界値問題の可視化

熱方程式の初期値境界値問題

$$u_t = u_{x x} \mbox{($(x,t)\in (0,1)\times (0,\infty)$)}$$

$$u(0,t) = u_(1,t)=0 \mbox{($t\in(0,\infty)$)}$$

$$u(x,0) = f(x) \mbox{($x\in[0,1]$)}$$

を差分法で解く。

/* heat1d.c -- 陰的スキーム (いわゆるθ法) で熱方程式を解く
*    空間1次元、同次 Dirichlet 境界条件の問題
*
*   「微分方程式と計算機演習」第11章 p237 のプログラムを修正・拡張したもの
C     ** IMPLICIT FINITE DIFFERENCE METHOD **
C     **         FOR HEAT EQUATION         **
C     nfunc : 関数の番号（5種類の関数を用意している）
C     theta : スキームのパラメーター
C     N : 分割の数
C     t : 時刻
C     tau : 時間の刻み幅
C     Tmax : 時刻の上限
C     h : 空間の刻み幅
C     AL,AD,AU : 係数行列
C     u : DIMENSION FOR UNKNOWN FUNCTION                     */

/*
*   このプログラムのコピーが欲しい場合は
*      cp /usr/local/meiji/D97/heat1d.c .
* とすればよい。
*
* 細かい拡張
*   1) 区間の左端、右端の座標を変数 a, b に設定するようにした。
*   2) グラフを描くごとに、それまでに描いたグラフを消去するかどうか、
*      変数 erase_always で制御するようにした。
*   3) 計算した数値データを表示するか変数 print_always で制御するよう
*      にした。
*   4) 何ステップおきにグラフを描き換えるか、数値データを表示するか、
*      変数 skips で制御するようにした。skips == 1 だと、毎回グラフを
*      書き換える。
*   5) 初期条件を与える関数 f(x, nfunc) で、登録されている関数の種類
*      を増やした。nfunc == 5 まである。特に nfunc == 4,5 は非対称な
*      関数である。
*
*/

#include <stdio.h>
#include <math.h>
#include <fplot.h>

#define ndim    1000

/* 円周率 (math.h にある M_PI の値を利用する
* あるいは double PI; として、どこかで PI = 4.0 * atan(1.0); と代入
*/

#define PI  M_PI

void print100(int, double, double *);
void axis(double, double, double, double);
void fsymbol(double, double, char *);
void print_data(double, double, double, double,
		int, double, int, double, double, double, double, double);
void save_picture(void);
void trid(int, double *, double *, double *, double *f);
void trilu(int n, double *al, double *ad, double *au);
void trisol(int n, double *al, double *ad, double *au, double *f);
double f(double, int);

int main()
{
    int N, nfunc, i, j, Jmax;
    double a, b, theta, h, Tmax, tau, c1, c2, c3, c4, lambda, t;
    double AL[ndim], AD[ndim], AU[ndim], ff[ndim], u[ndim + 1];
    double xleft, ybottom, xright, ytop;
/* 何ステップおきにグラフを書き換えるか */
    int skips = 1;
    double dt;
/* 以前描いたグラフを消すか？(0=No, 1=Yes) */
    int erase_always = 0;
/* 数値を表示するか?(0=No, 1=Yes) */
    int print_always = 0;

/* 区間の左端、右端 */
    a = 0.0;
    b = 1.0;

/* 入力 */
    printf("入力して下さい : nfunc(1..5)=");
    scanf("%d", &nfunc);
    printf("入力して下さい : θ=");
    scanf("%lf", &theta);
    printf("入力して下さい : N(<=%d)=", ndim);
    scanf("%d", &N);
    if (N <= 1 || N > ndim)
	return 0;
    printf("入力して下さい : λ=");
    scanf("%lf", &lambda);

    h = (b - a) / N;
    tau = lambda * h * h;
    printf(" 時間の刻み幅τ= %g になりました。\n", tau);

    printf("入力して下さい : Tmax=");
    scanf("%lf", &Tmax);

    printf("図を描く時間間隔 Δt は : ");
    scanf("%lf", &dt);
    skips = (dt + 0.1 * tau) / tau;

/* 係数行列の計算 */
    c1 = -theta * lambda;
    c2 = 1.0 + 2.0 * theta * lambda;
    c3 = 1.0 - 2.0 * (1.0 - theta) * lambda;
    c4 = (1.0 - theta) * lambda;
    for (i = 1; i < N; i++) {
	AL[i] = c1;
	AD[i] = c2;
	AU[i] = c1;
    }

/* LU 分解する */
    trilu(N - 1, AL + 1, AD + 1, AU + 1);

/* 初期値 */
    for (i = 0; i <= N; i++)
	u[i] = f(a + i * h, nfunc);

/* 出力 (t=0) */
/* 出力 (ここでは画面に数値を表示すること) は Fortran と C で
   かなり異なるので、直接の翻訳は出来ない。ごたごたするので、
   print100 という関数にまとめた */
    t = 0.0;
    print100(N, t, u);

/* グラフィックス画面の準備 -- まず画面の範囲を決めてから */
    xleft = a - (b - a) / 10;
    xright = b + (b - a) / 10;
    ybottom = -0.1;
    ytop = 1.1;

/* fplot ライブラリィの呼び出し */
    openpl();
    fspace2(xleft, ybottom, xright, ytop);

/* パラメーター、入力データ等の表示 */
    print_data(xleft, ybottom, xright, ytop,
	       nfunc, theta, N, lambda, tau, Tmax, t, dt);

/* 解 u の t=0 におけるグラフを描く */
    fmove(a, u[0]);
    for (i = 1; i <= N; i++)
	fcont(a + i * h, u[i]);
    xsync();

    Jmax = (Tmax + 0.1 * tau) / tau;

/* 繰り返し計算 */
    for (j = 1; j <= Jmax; j++) {

	/* 連立一次方程式を作って解く */
	for (i = 1; i < N; i++)
	    ff[i] = c3 * u[i] + c4 * (u[i - 1] + u[i + 1]);
	trisol(N - 1, AL + 1, AD + 1, AU + 1, ff + 1);

	/* 求めた値を u[] に収める */
	for (i = 1; i < N; i++)
	    u[i] = ff[i];
	u[0] = 0.0;
	u[N] = 0.0;

	/* 時刻の計算 */
	t = j * tau;

	if (j % skips == 0) {
	    /* 数値データの表示 */
	    if (print_always)
		print100(N, t, u);

	    /* 解のグラフを描く */
	    if (erase_always) {
		/* これまで描いたものを消し、パラメーターを再表示する */
		erase();
		print_data(xleft, ybottom, xright, ytop,
			   nfunc, theta, N, lambda, tau, Tmax, t, dt);
	    }
	    fmove(a, u[0]);
	    for (i = 1; i <= N; i++)
		fcont(a + i * h, u[i]);
	    xsync();
	}
    }

    save_picture();

    printf("終了したければ、ウィンドウをマウスでクリックして下さい。\n");
    closepl();
    return 0;
}

/**********************************************************/

/* 初期条件を与える関数。複数登録して番号 nfunc で選択する。 */

double f(double x, int nfunc)
{
/* f(x)=max(x,1-x) */
    if (nfunc == 1) {
	if (x <= 0.5)
	    return x;
	else
	    return 1.0 - x;
    }
/* f(x)=1 */
    else if (nfunc == 2)
	return 1.0;
    else if (nfunc == 3)
	return sin(PI * x);
    /* f(x)= 変な形をしたもの(by M.Sakaue) */
    else if (nfunc == 4) {
	if (x <= 0.1)
	    return 5 * x;
	else if (x <= 0.3)
	    return -2 * x + 0.7;
	else if (x <= 0.5)
	    return 4.5 * x - 1.25;
	else if (x <= 0.7)
	    return -x + 1.5;
	else if (x <= 0.9)
	    return x + 0.1;
	else
	    return -10 * x + 10.0;
    }
/* やはり非対称な形をしたもの(by M.Sakaue) */
    else if (nfunc == 5) {
	if (x <= 0.2)
	    return -x + 0.8;
	else if (x <= 0.3)
	    return 4 * x - 0.2;
	else if (x <= 0.4)
	    return -4 * x + 2.2;
	else if (x <= 0.7)
	    return -x + 1.0;
	else if (x <= 0.8)
	    return 4 * x - 2.6;
	else
	    return -x + 1.5;
    }
}

/**********************************************************/

/* 配列 u[] の内容(u[0],u[1],...,u[n])を一行 5 個ずつ表示する。 */

void print100(int N, double t, double *u)
{
    int i;
    printf("T= %12.4e\n", t);
    for (i = 0; i < 5; i++)
	printf("  I     u(i)    ");
    printf("\n");
    for (i = 0; i <= N; i++) {
	printf("%3d%12.4e ", i, u[i]);
	if (i % 5 == 4)
	    printf("\n");
    }
    printf("\n");
}

/**********************************************************/

/* ウィンドウに座標軸を表示する */
void axis(double xleft, double ybottom, double xright, double ytop)
{
    linemod("dotted");
    fline(xleft, 0.0, xright, 0.0);
    fline(0.0, ybottom, 0.0, ytop);
    linemod("solid");
}

/* ウィンドウ内の位置 (x,y) から文字列 s を表示する */
void fsymbol(double x, double y, char *s)
{
    fmove(x, y);
    label(s);
}

#define X(x)    (xleft+(x)*(xright-xleft))
#define Y(y)    (ybottom+(y)*(ytop-ybottom))

/* パラメーターの値等をウィンドウに表示する */
void print_data(double xleft, double ybottom, double xright, double ytop,
		int nfunc, double theta, int N, double lambda, double tau,
		double Tmax, double t, double dt)
{
    char message[80];

    axis(xleft, ybottom, xright, ytop);
    sprintf(message, "heat equation, u(0)=u(1)=0");
    fsymbol(X(0.2), Y(0.95), message);
    sprintf(message,
	    "nfunc=%d, theta=%g, N=%d, lambda=%g, tau=%g, Tmax=%g, dt=%g",
	    nfunc, theta, N, lambda, tau, Tmax, dt);
    fsymbol(X(0.2), Y(0.9), message);
    sprintf(message, "t = %g", t);
    if (t != 0.0)
	fsymbol(X(0.2), Y(0.85), message);
}

/* 現在ウィンドウに表示されている図をファイルに保存する */
void save_picture()
{
    char filename[200];
    printf("図を保存するファイル名: ");
    scanf("%s", filename);
    mkplot(filename);
}

/*************************************************************
 *************************************************************
 *************************************************************/

/*
 *  三重対角行列に対する LU 分解に基づく連立1次方程式の求解
 *
 *  (Gauss の消去法を用いているが、ピボットの選択等はしていないので、
 *   係数行列が正定値対称行列でないと結果は保証されない。)
 *
 *  n は未知数の個数。
 *  f は最初は連立1次方程式の右辺のベクトル b。最後は解ベクトル x。
 *   (添字は 0 から。i.e. b[0],b[1],...,b[n-1] にデータが入っている。)
 *
 *  AL, AD, AU は係数行列の
 *     下三角部分 (lower part)
 *     対角部分 (diagonal part)
 *     上三角部分 (upper part)
 *  つまり
 *
 *  AD[0] AU[0]   0   .................. 0
 *  AL[1] AD[1] AU[1]   0  ............. 0
 *    0   AL[2] AD[2] AU[2]  0 ......... 0
 *          ....................
 *                    AL[n-2] AD[n-2] AU[n-2]
 *                      0     AL[n-1] AD[n-1]
 *
 *  ここで AL[0], AU[n-1] は意味がないことに注意。
 */

/* 三項方程式(係数行列が三重対角である連立一次方程式のこと)を解く
 *   入力
 *     n: 未知数の個数
 *     al,ad,au: 連立一次方程式の係数行列
 *       (al: 対角線の下側、 ad: 対角線、 au: 対角線の上側)
 *        al[i] = A_{i,i-1}, ad[i] = A_{i,i}, au[i] = A_{i,i+1},
 *        al[0], au[n-1] は意味がない)
 *        f: 連立一次方程式の右辺の既知ベクトル b
 *   出力
 *     al,ad,au: 入力した係数行列を LU 分解したもの
 *     f: 連立一次方程式の解 x
 *   能書き
 *     一度 call すると係数行列を LU 分解したものが返されるので、
 *     以後は同じ係数行列に関する連立一次方程式を解くために、
 *     サブルーチン trisol が使える。
 *   注意
 *     ピボットの選択をしていないので、係数行列が正定値である
 *     などの適切な条件がない場合は結果が保証できない。
 */

/* 三項方程式を解く */
void trid(int n, double *al, double *ad, double *au, double *f)
{
    trilu(n, al, ad, au);
    trisol(n, al, ad, au, f);
}

/* 三重対角行列の LU 分解 (pivoting なし) */
void trilu(int n, double *al, double *ad, double *au)
{
    int i, nm1 = n - 1;
    /* 前進消去 (forward elimination) */
    for (i = 0; i < nm1; i++)
	ad[i + 1] -= au[i] * al[i + 1] / ad[i];
}

/* LU 分解済みの三重対角行列を係数に持つ三項方程式を解く */
void trisol(int n, double *al, double *ad, double *au, double *f)
{
    int i, nm1 = n - 1;
    /* 前進消去 (forward elimination) */
    for (i = 0; i < nm1; i++)
	f[i + 1] -= f[i] * al[i + 1] / ad[i];
    /* 後退代入 (backward substitution) */
    f[nm1] /= ad[nm1];
    for (i = n - 2; i >= 0; i--)
	f[i] = (f[i] - au[i] * f[i + 1]) / ad[i];
}