5.3 陽解法 heat1d-e-eggx.c

/*
 * heat1d-e-eggx.c -- 1次元熱伝導方程式の初期値境界値問題を陽解法で解く。
 *   コンパイル:  egg heat1d-e-eggx.c -o heat1d-e-eggx
 *
 * オリジナルは fplot ライブラリィを利用した
 *  http://nalab.mind.meiji.ac.jp/%7Emk/program/fdm/heat1d-e.c
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <eggxlib.h>

int main()
{
  int i, n, nMax, N;
  double tau, h, lambda, Tmax;
  double *u, *newu;
  double f(double);
  int win;
  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));
  
  /* 初期値の代入 */
  for (i = 0; i <= N; i++)
    u[i] = f(i * h);

  /* ***************** グラフィックスの準備 ***************** */
  win = eggx_gopen(600, 600);
  eggx_gsetbgcolor(win, "black"); eggx_gclr(win);
  eggx_newcolor(win, "green");
  /* 座標系の定義: [-0.1,1.1]×[-0.1,1.1] という閉領域を表示する */
  eggx_window(win, -0.1, -0.1, 1.1, 1.1);

  /* タイトルと入力パラメーターを表示する */
  eggx_drawstr(win, 0.1, 0.8, 14, 0,
         "heat equation, homogeneous Dirichlet boundary condition");
  sprintf(message, "N=%d, lambda=%g, Tmax=%g", N, lambda, Tmax);
  eggx_drawstr(win, 0.1, 0.7, 14, 0, message);

  /* 座標軸を表示する */
  eggx_line(win, -0.1, 0.0, PENUP);
  eggx_line(win,  1.1, 0.0, PENDOWN);
  eggx_line(win,  0.0, -0.1, PENUP);
  eggx_line(win,  0.0, 1.1, PENDOWN);

  /* t=0 の状態を表示する */
  eggx_line(win, 0.0, u[0], PENUP);
  for (i = 1; i <= N; i++)
    eggx_line(win, i * h, u[i], PENDOWN);

  /* ループを何回まわるか計算する (四捨五入) */
  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)τ) の状態を表示する */
    eggx_line(win, 0.0, u[0], PENUP);
    for (i = 1; i <= N; i++)
      eggx_line(win, i * h, u[i], PENDOWN);
  }
  
  printf("終りました。ウィンドウ内で文字を入力してください。\n");
  eggx_ggetch(win);
  /* ウィンドウを閉じる */
  eggx_gclose(win);
  return 0;
}

double f(double x)
{
  if (x <= 0.5)
    return x;
  else
    return 1.0 - x;
}