10.1 陽解法

// heat1d-e-freefem.edp --- 熱方程式を差分法 (前進Euler法) で解く
// 実行: 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);
}