1 /*
2 * heat1d-i-glsc.c -- 陰的スキームで熱方程式を解く
3 *
4 * 菊地・山本, 微分方程式と計算機演習, 山海堂 (1991) 第11章
5 * p237 のプログラムを修正・拡張したもの
6 *
7 C ** IMPLICIT FINITE DIFFERENCE METHOD **
8 C ** FOR HEAT EQUATION **
9 C nfunc : 関数の番号（5種類の関数を用意している）
10 C theta : スキームのパラメーター
11 C N : 分割の数
12 C t : 時刻
13 C tau : 時間の刻み幅
14 C Tmax : 時刻の上限
15 C h : 空間の刻み幅
16 C AL,AD,AU : 係数行列
17 C u : DIMENSION FOR UNKNOWN FUNCTION
18 *
19 * このプログラムは
20 * http://www.math.meiji.ac.jp/%7Emk/program/
21 * から入手可能
22 */
23
24 #include <stdio.h>
25 #include <stdlib.h>
26 #include <math.h>
27
28 #define G_DOUBLE
29 #include <glsc.h>
30
31 /* 円周率 --- main() で pi = 4.0 * atan(1.0); と代入 */
32
33 double pi;
34
35 void trilu(int, double *, double *, double *);
36 void trisol(int, double *, double *, double *, double *);
37 double f(double, int);
38 void print100(int, double, double *);
39 void axis(double, double, double, double);
40 void fsymbol(double, double, char *);
41 void print_parameters(double, double, double, double,
42 int, double, int,
43 double, double, double, double);
44 void save_picture(void);
45
46 int main()
47 {
48 int N, nfunc, i, n, nMax;
49 double a, b, theta, h, Tmax, tau, c1, c2, c3, c4, lambda, t;
50 double *AL, *AD, *AU, *ff, *u;
51 double xleft, ybottom, xright, ytop;
52 /* グラフを書き換える時間間隔 */
53 double dt;
54 /* 何ステップおきにグラフを書き換えるか */
55 int skip;
56 /* 以前描いたグラフを消すか？ (0=No, 1=Yes) */
57 int erase_always = 0;
58 /* 数値を表示するか? (0=No, 1=Yes) */
59 int print_numerical_data = 0;
60 /* */
61 double win_width, win_height, w_margin, h_margin;
62
63 /* 円周率 */
64 pi = 4.0 * atan(1.0);
65 /* 区間の左端、右端 */
66 a = 0.0;
67 b = 1.0;
68
69 /* 入力 */
70 printf("入力して下さい : nfunc(1..5)=");
71 scanf("%d", &nfunc);
72 printf("入力して下さい : θ=");
73 scanf("%lf", &theta);
74 printf("入力して下さい : N=");
75 scanf("%d", &N);
76 printf("入力して下さい : λ=");
77 scanf("%lf", &lambda);
78
79 /* ベクトルの準備 */
80 if ((AL = malloc(N * sizeof(double))) == NULL) {
81 fprintf(stderr, "cannot allocate AL[]\n");
82 exit(1);
83 }
84 if ((AD = malloc(N * sizeof(double))) == NULL) {
85 fprintf(stderr, "cannot allocate AD[]\n");
86 exit(1);
87 }
88 if ((AU = malloc(N * sizeof(double))) == NULL) {
89 fprintf(stderr, "cannot allocate AU[]\n");
90 exit(1);
91 }
92 if ((ff = malloc(N * sizeof(double))) == NULL) {
93 fprintf(stderr, "cannot allocate ff[]\n");
94 exit(1);
95 }
96 if ((u = malloc((N + 1) * sizeof(double))) == NULL) {
97 fprintf(stderr, "cannot allocate u[]\n");
98 exit(1);
99 }
100 /* */
101 h = (b - a) / N;
102 tau = lambda * h * h;
103 printf(" 時間の刻み幅τ= %g になりました。\n", tau);
104
105 /* 入力 */
106 printf("入力して下さい : 最終時刻 Tmax=");
107 scanf("%lf", &Tmax);
108 printf("入力して下さい : グラフ書き換え時間間隔 Δt=");
109 scanf("%lf", &dt);
110 skip = (int) rint(dt / tau);
111 /* skip が 0 になっていないかチェック */
112 if (skip == 0) {
113 fprintf(stderr, "Δtが小さすぎるので、Δt=τ=%g にします。\n", tau);
114 skip = 1;
115 dt = tau;
116 }
117
118 /* 係数行列の準備 */
119 c1 = - theta * lambda;
120 c2 = 1.0 + 2.0 * theta * lambda;
121 c3 = 1.0 - 2.0 * (1.0 - theta) * lambda;
122 c4 = (1.0 - theta) * lambda;
123 for (i = 1; i < N; i++) {
124 AL[i] = c1;
125 AD[i] = c2;
126 AU[i] = c1;
127 }
128
129 /* 初期値 */
130 for (i = 0; i <= N; i++)
131 u[i] = f(a + i * h, nfunc);
132
133 /* 出力（t=0） */
134 /* 出力（ここでは画面に数値を表示すること）は FORTRAN と C で
135 かなり異なるので、直接の翻訳は出来ない。ごたごたするので、
136 print100 という関数にまとめた */
137 t = 0.0;
138 print100(N, t, u);
139
140 /* ***************** グラフィックスの準備 ***************** */
141 /* グラフィックス画面の準備 -- まず画面の範囲を決めてから */
142 xleft = a - (b - a) / 10;
143 xright = b + (b - a) / 10;
144 ybottom = - 0.1;
145 ytop = 1.1;
146 /* メタファイル名は "HEAT",
147 * ウィンドウのサイズは、
148 横 win_width + 2 * w_margin, 縦 win_height + 2 * h_margin */
149 win_width = 200.0; win_height = 200.0; w_margin = 10.0; h_margin = 10.0;
150 g_init("HEAT", win_width + 2 * w_margin, win_height + 2 * h_margin);
151 /* 画面とメタファイルの両方に記録する */
152 g_device(G_BOTH);
153 /* 座標系の定義: [-0.1,1.1]×[-0.1,1.1] という閉領域を表示する */
154 g_def_scale(0,
155 -0.1, 1.1, -0.1, 1.1,
156 w_margin, h_margin, win_width, win_height);
157 /* 線を二種類用意する */
158 g_def_line(0, G_BLACK, 0, G_LINE_SOLID);
159 g_def_line(1, G_BLACK, 0, G_LINE_DOTS);
160 /* 表示するための文字列の属性を定義する */
161 g_def_text(0, G_BLACK, 3);
162 /* 定義したものを選択する */
163 g_sel_scale(0); g_sel_line(0); g_sel_text(0);
164
165 /* タイトルと入力パラメーターを表示する */
166 print_parameters(xleft, ybottom, xright, ytop,
167 nfunc, theta, N,
168 lambda, tau, Tmax, t);
169
170 /* 座標軸を表示する */
171 g_sel_line(1);
172 g_move(-0.1, 0.0); g_plot(1.1, 0.0);
173 g_move(0.0, -0.1); g_plot(0.0, 1.1);
174 g_sel_line(0);
175
176 /* t=0 の状態を表示する */
177 g_move(0.0, u[0]);
178 for (i = 1; i <= N; i++)
179 g_plot(i * h, u[i]);
180
181 /* 係数行列の LU 分解 */
182 trilu(N - 1, AL + 1, AD + 1, AU + 1);
183
184 /* 繰り返し計算 */
185 nMax = (int) rint(Tmax / tau);
186 for (n = 1; n <= nMax; n++) {
187 /* 連立1次方程式を作って解く */
188 for (i = 1; i < N; i++)
189 ff[i] = c3 * u[i] + c4 * (u[i - 1] + u[i + 1]);
190 trisol(N - 1, AL + 1, AD + 1, AU + 1, ff + 1);
191 for (i = 1; i < N; i++)
192 u[i] = ff[i];
193 /* 境界条件 */
194 u[0] = 0.0;
195 u[N] = 0.0;
196
197 /* 時刻 */
198 t = n * tau;
199
200 if (n % skip == 0) {
201 /* 数値データの表示 */
202 if (print_numerical_data)
203 print100(N, t, u);
204
205 /* 解のグラフを描く */
206 if (erase_always) {
207 /* これまで描いたものを消し、パラメーターを再表示する */
208 g_cls();
209 print_parameters(xleft, ybottom, xright, ytop,
210 nfunc, theta, N, lambda, tau, Tmax, t);
211 }
212 g_move(0.0, u[0]);
213 for (i = 1; i <= N; i++)
214 g_plot(i * h, u[i]);
215 }
216 }
217
218 printf("マウスでウィンドウをクリックして下さい。\n");
219 g_sleep(-1.0);
220 g_term();
221
222 return 0;
223 }
224
225 /**********************************************************/
226
227 /* 三重対角行列の LU 分解 (pivoting なし) */
228 void trilu(int n, double *al, double *ad, double *au)
229 {
230 int i, nm1 = n - 1;
231 /* 前進消去 (forward elimination) */
232 for (i = 0; i < nm1; i++) {
233 al[i + 1] /= ad[i];
234 ad[i + 1] -= au[i] * al[i + 1];
235 }
236 }
237
238 /* LU 分解済みの三重対角行列を係数に持つ三項方程式を解く */
239 void trisol(int n, double *al, double *ad, double *au, double *b)
240 {
241 int i, nm1 = n - 1;
242 /* 前進消去 (forward elimination) */
243 for (i = 0; i < nm1; i++) b[i + 1] -= b[i] * al[i + 1];
244 /* 後退代入 (backward substitution) */
245 b[nm1] /= ad[nm1];
246 for (i = n - 2; i >= 0; i--) b[i] = (b[i] - au[i] * b[i + 1]) / ad[i];
247 }
248
249 /**********************************************************/
250
251 /* 初期条件を与える関数。複数登録して番号 nfunc で選択する。 */
252
253 double f(double x, int nfunc)
254 {
255 /* f(x)=max(x,1-x) */
256 if (nfunc == 1) {
257 if (x <= 0.5)
258 return x;
259 else
260 return 1.0 - x;
261 }
262 /* f(x)=1 */
263 else if (nfunc == 2)
264 return 1.0;
265 else if (nfunc == 3)
266 return sin(pi * x);
267 /* f(x)= 変な形をしたもの(by M.Sakaue) */
268 else if (nfunc == 4) {
269 if (x <= 0.1)
270 return 5 * x;
271 else if (x <= 0.3)
272 return -2 * x + 0.7;
273 else if (x <= 0.5)
274 return 4.5 * x - 1.25;
275 else if (x <= 0.7)
276 return -x + 1.5;
277 else if (x <= 0.9)
278 return x + 0.1;
279 else
280 return -10 * x + 10.0;
281 }
282 /* やはり非対称な形をしたもの(by M.Sakaue) */
283 else if (nfunc == 5) {
284 if (x <= 0.2)
285 return -x + 0.8;
286 else if (x <= 0.3)
287 return 4 * x - 0.2;
288 else if (x <= 0.4)
289 return -4 * x + 2.2;
290 else if (x <= 0.7)
291 return -x + 1.0;
292 else if (x <= 0.8)
293 return 4 * x - 2.6;
294 else
295 return -x + 1.5;
296 } else
297 return 0;
298 }
299
300 /**********************************************************/
301
302 /* 配列 u[] の内容(u[0],u[1],...,u[n])を一行 5 個ずつ表示する。 */
303
304 void print100(int N, double t, double *u)
305 {
306 int i;
307 printf("T= %12.4e\n", t);
308 for (i = 0; i < 5; i++)
309 printf(" I u(i) ");
310 printf("\n");
311 for (i = 0; i <= N; i++) {
312 printf("%3d%12.4e ", i, u[i]);
313 if (i % 5 == 4)
314 printf("\n");
315 }
316 printf("\n");
317 }
318
319 /**********************************************************/
320
321 /* ウィンドウに座標軸を表示する */
322 void axis(double xleft, double ybottom, double xright, double ytop)
323 {
324 g_sel_line(1);
325 g_move(xleft, 0.0); g_plot(xright, 0.0);
326 g_move(0.0, ybottom); g_plot(0.0, ytop);
327 g_sel_line(0);
328 }
329
330 /* パラメーターの値等をウィンドウに表示する */
331 void print_parameters(double xleft, double ybottom, double xright, double ytop,
332 int nfunc, double theta, int N,
333 double lambda, double tau, double Tmax, double t)
334 {
335 char message[80];
336
337 axis(xleft, ybottom, xright, ytop);
338
339 g_text(30.0, 30.0, "heat equation, u(0)=u(1)=0");
340 sprintf(message, "nfunc=%d, theta=%g, N=%d, lambda=%g, tau=%g, Tmax=%g",
341 nfunc, theta, N, lambda, tau, Tmax);
342 g_text(30.0, 60.0, message);
343 if (t != 0.0) {
344 sprintf(message, "t = %g", t);
345 g_text(30.0, 90.0, message);
346 }
347
348 }