/* * heat1d-e-glsc.c -- 1次元熱伝導方程式の初期値境界値問題 * * u_t(x,t)=u_{xx}(x,t) x∈(0,1), t∈(0,∞) * u(0,t)=u(1,t)=0 t∈(0,∞) * u(x,0)=f(x) x∈[0,1] * * を陽解法で解く。 * * コンパイル: ccmg heat1d-e-glsc.c * 文字コードに注意！ Windows では qkc -ms heat1d-e-glsc.c が必要かも * * オリジナルは fplot ライブラリィを利用した * http://nalab.mind.meiji.ac.jp/~mk/program/fdm/heat1d-e.c */ #include #include #include #ifndef G_DOUBLE #define G_DOUBLE #endif #include int main() { int i, n, nMax, N; double tau, h, lambda, Tmax; double *u, *newu; double f(double); double win_width, win_height, w_margin, h_margin; char message[100]; /* N, λ を入力する */ printf("区間の分割数 N = "); scanf("%d", &N); printf("λ (=τ/h^2) = "); scanf("%lf", &lambda); /* 時間の刻み幅 h, τ を計算する */ h = 1.0 / N; tau = lambda * h * h; printf("τ=%g\n", tau); /* 最終時刻を入力する */ printf("最終時刻 Tmax = "); scanf("%lf", &Tmax); /* ベクトル u, newu を用意する */ u = malloc(sizeof(double) * (N+1)); newu = malloc(sizeof(double) * (N+1)); if (u == NULL || newu == NULL) { fprintf(stderr, "メモリーが準備できませんでした。\n"); exit(1); } /* 初期値の代入 */ for (i = 0; i <= N; i++) u[i] = f(i * h); /* ***************** グラフィックスの準備 ***************** */ /* メタファイル名は "HEAT1DE", * ウィンドウのサイズは、 横 win_width + 2 * w_margin, 縦 win_height + 2 * h_margin */ win_width = 150.0; win_height = 150.0; w_margin = 10.0; h_margin = 10.0; g_init("HEAT1DE", win_width + 2 * w_margin, win_height + 2 * h_margin); /* 画面とメタファイルの両方に記録する */ g_device(G_BOTH); /* 座標系の定義: [-0.1,1.1]×[-0.1,1.1] という閉領域を表示する */ g_def_scale(0, -0.1, 1.1, -0.1, 1.1, w_margin, h_margin, win_width, win_height); /* 線を二種類用意する */ g_def_line(0, G_BLACK, 0, G_LINE_SOLID); g_def_line(1, G_BLACK, 0, G_LINE_DOTS); /* 表示するための文字列の属性を定義する */ g_def_text(0, G_BLACK, 3); /* 定義したものを選択する */ g_sel_scale(0); g_sel_line(0); g_sel_text(0); /* タイトルと入力パラメーターを表示する */ g_text(30.0, 30.0, "heat equation, homogeneous Dirichlet boundary condition"); snprintf(message, sizeof(message), "N=%d, lambda=%g, Tmax=%g", N, lambda, Tmax); g_text(30.0, 60.0, message); /* 座標軸を表示する */ g_sel_line(1); g_move(-0.1, 0.0); g_plot(1.1, 0.0); g_move(0.0, -0.1); g_plot(0.0, 1.1); g_sel_line(0); /* t=0 の状態を表示する */ g_move(0.0, u[0]); for (i = 1; i <= N; i++) g_plot(i * h, u[i]); /* ループを何回まわるか計算する (四捨五入) */ nMax = rint(Tmax / tau); /* 時間に関するステップを進めるループ */ for (n = 0; n < nMax; n++) { /* 差分方程式 (n -> n+1) */ for (i = 1; i < N; i++) newu[i] = (1.0 - 2 * lambda) * u[i] + lambda * (u[i+1] + u[i-1]); /* 計算値を更新 */ for (i = 1; i < N; i++) u[i] = newu[i]; /* Dirichlet 境界条件 */ u[0] = u[N] = 0.0; /* この時刻 (t=(n+1)τ) の状態を表示する */ g_move(0.0, u[0]); for (i = 1; i <= N; i++) g_plot(i * h, u[i]); } printf("終りました。X の場合はウィンドウをクリックして下さい。\n"); g_sleep(-1.0); /* ウィンドウを閉じる */ g_term(); return 0; } /* 初期条件 --- ここでは min(x,1-x) という山形のグラフを持つ関数 */ double f(double x) { if (x <= 0.5) return x; else return 1.0 - x; }