/* 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 #include #define ndim 1000 /* 円周率 (math.h にある M_PI の値を利用する * あるいは double PI; として、どこかで PI = 4.0 * atan(1.0); と代入 */ #define PI M_PI main() { int N, nfunc, i, j, Jmax; double a, b, theta, h, Tmax, tau, c1, c2, c3, c4, lambda, f(), 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(x, nfunc) 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 個ずつ表示する。 */ print100(N, t, u) int N; double t, 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"); } /**********************************************************/ /* ウィンドウに座標軸を表示する */ axis(xleft, ybottom, xright, ytop) double xleft, ybottom, xright, ytop; { linemod("dotted"); fline(xleft, 0.0, xright, 0.0); fline(0.0, ybottom, 0.0, ytop); linemod("solid"); } /* ウィンドウ内の位置 (x,y) から文字列 s を表示する */ fsymbol(x, y, s) double x, y; char s[]; { fmove(x, y); label(s); } #define X(x) (xleft+(x)*(xright-xleft)) #define Y(y) (ybottom+(y)*(ytop-ybottom)) /* パラメーターの値等をウィンドウに表示する */ print_data(xleft, ybottom, xright, ytop, nfunc, theta, N, lambda, tau, Tmax, t, dt) double xleft, ybottom, xright, ytop, theta, lambda, tau, Tmax, t, dt; int nfunc, N; { 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); } /* 現在ウィンドウに表示されている図をファイルに保存する */ 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 が使える。 * 注意 * ピボットの選択をしていないので、係数行列が正定値である * などの適切な条件がない場合は結果が保証できない。 */ /* 三項方程式を解く */ trid(n,al,ad,au,f) int n; double al[],ad[],au[],f[]; { trilu(n,al,ad,au); trisol(n,al,ad,au,f); } /* 三重対角行列の LU 分解 (pivoting なし) */ trilu(n,al,ad,au) int n; double al[],ad[],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 分解済みの三重対角行列を係数に持つ三項方程式を解く */ trisol(n,al,ad,au,f) int n; double al[],ad[],au[],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]; }