B..1 例: 差分法で解いた1次元熱伝導方程式の初期値境界値問題の解の可視化

1次元区間 $(0,1)=\left\{x\in\R\relmiddle\vert 0<x<1\right\}$ における熱伝導方程式の初期値境界値問題

$$\frac{\rd u}{\rd t}(x,t) =\frac{\rd^2 u}{\rd x^2}(x,t) x\in (0,1) t\in(0,T_{\mathrm{max}}) ,$$

$$u(0,t)=u(1,t)=0 t\in(0,T_{\mathrm{max}}) ,$$

$$u(x,0)=u_0(x) x\in[0,1] .$$

を差分法 (陽解法) で解いて、それを可視化する。

ここで $u_0$ は時刻 $t=0$ での初期温度分布である。

時間発展する系なので、 各時刻 $t=t_n$ での差分解 $x\mapsto u(x,t_n)$ のグラフを描くことにする。

差分法の説明は、例えば 「発展系の数値解析」 に載っている。

/*
 * heat1d-gnuplot.c --- 1次元領域の熱伝導方程式
 *   u_t=u_{xx}, 同次Dirichlet境界条件
 */

#include <stdio.h>
#include <math.h>   // floor()
#include <stdlib.h> // malloc()
#include <unistd.h> // sleep()

// 初期条件
double u0(double x)
{
  if (x < 0.5)
    return x;
  else
    return 1.0 - x;
}

int main(void)
{
  int N, i, n, nMax, nskip;
  double Tmax, tau, h, t, lambda, c, dt;
  double *x, *u, *newu;
  FILE *gp;
  Tmax = 0.5; // 最終時刻
  dt = 0.01; // 描画間隔
  // 分割
  N = 50;
  h = 1.0 / N;
  lambda = 0.5; // λ=τ/h^2
  tau = lambda * h * h;
  nskip = rint(dt / tau); if (nskip == 0) nskip = 1;
  nMax = Tmax / tau;
  printf("N=%d, h=%g, tau=%g, nskip=%d\n", N, h, tau, nskip);
  // メモリー確保
  x = malloc(sizeof(double) * (N+1));
  u = malloc(sizeof(double) * (N+1));
  newu = malloc(sizeof(double) * (N+1));
  if (x == NULL || u == NULL || newu == NULL) {
    fprintf(stderr, "メモリー確保に失敗しました。\n");
    exit(1);
  }
  // gnuplot を呼び出す
  if ((gp = popen("gnuplot -persist", "w")) == NULL) {
    fprintf(stderr, "gnuplot 呼び出しに失敗しました。\n");
    exit(1);
  }
  // fprintf(gp, "set term gif animate\n");
  // fprintf(gp, "set output \"heat1d.gif\"\n");
  fprintf(gp, "set yrange [-0.01:0.5]\n");
  fprintf(gp, "set xlabel \"x\"\n");
  fprintf(gp, "set ylabel \"u\"\n");
  fprintf(gp, "set label \"1D heat equation\" at screen 0.4,0.95\n");
  // 初期値
  t = 0;
  for (i = 0; i <= N; i++) {
    x[i] = i * h;
    u[i] = u0(x[i]);
  }
  // 初期値を描く
  fprintf(gp, "plot '-' with lines title \"t=0\"\n");
  for (i = 0; i <= N; i++) {
      fprintf(gp, "%f %f\n", x[i], u[i]);
  }
  fprintf(gp, "e\n"); fflush(gp);
  sleep(3);
  // 時間を進める
  c = 1 - 2.0 * lambda;
  for (n = 1; n <= nMax; n++) {
    t = n * tau;
    for (i = 1; i < N; i++)
        newu[i] = c * u[i] + lambda * (u[i+1]+u[i-1]);
    newu[0] = newu[N] = 0.0;
    for (i = 0; i <= N; i++)
        u[i] = newu[i];
    if (n % nskip == 0) {
      fprintf(gp, "plot '-' with lines title \"t=%g\"\n", t);
      for (i = 0; i <= N; i++)
          fprintf(gp, "%f %f\n", x[i], u[i]);
      fprintf(gp, "e\n"); fflush(gp);
    }
  }
  pclose(gp);
  return 0;
}

Mac のターミナルの場合: プログラムの入手、コンパイル、実行

curl -O http://nalab.mind.meiji.ac.jp/~mk/misc/20191229/heat1d-gnuplot.c
cc heat1vi d-gnuplot.c
./a.out

(Xcode の cc, GNU の gcc で動作確認済み)

GIFアニメーション