heat1d-e-freefem.edp -- 1次元熱方程式に対する差分法プログラム |
// heat1d-e-freefem.edp // 実行: FreeFem++ heat1d-e-freefem.edp // N, x が大域的な識別子を隠すと警告が出る (つまり名前が衝突する)。 int i, N=50, n, nMax; real h = 1.0 / N, lambda = 0.5, Tmax = 1.0; real tau = lambda * h^2; real[int] x(N+1); real[int] u(N+1); real[int] newu(N+1); cout << "h= " << h << endl; cout << "lambda= " << lambda << endl; cout << "tau= " << tau << endl; func real f(real x) { if (x < 0.5) return x; else return 1 - x; } for (i = 0; i <= N; i++) { x[i] = i * h; u[i] = f(x[i]); } plot([x,u], bb=[[-0.1,-0.1],[1.1,1.1]],aspectratio=true, wait=true); nMax = rint(Tmax / tau); cout << "nMax =" << nMax << endl; for (n = 1; n <= nMax; n++) { for (i = 1; i < N; i++) newu[i] = (1.0 - 2.0 * lambda) * u[i] + lambda * (u[i+1] + u[i-1]); // cout << newu << endl; newu[0] = newu[N] = 0; u = newu; plot([x, u], bb=[[-0.1,-0.1],[1.1,1.1]],aspectratio=true); } |
なお、C++ ライクな cin も使うことが出来るので、 N の値をキーボード入力することも可能である。
int i, N, n, nMax; cout << "input N: "; cin >> N; real h = 1.0 / N, lambda = 0.5, Tmax = 1.0; |