7.3 陽解法 heat1d-e-glut.c

/*
 * heat1d-e-glut.c -- 1次元熱伝導方程式の初期値境界値問題を陽解法で解く。
 *   gcc -I/usr/local/include -o heat1d heat1d-e-glut.c -lglut -lGL -lX11 -lm
 *
 * オリジナルは fplot ライブラリィを利用した
 *  http://nalab.mind.meiji.ac.jp/%7Emk/program/fdm/heat1d-e.c
 */

#define WIDTH (600)
#define HEIGHT (600)

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <GL/glut.h>

double xc = 0.0, yc = 0.0, ratio_x = 0.9, ratio_y = 0.9;

void screen(double xl, double yb, double xr, double yt)
{
  /* [xl,xr]×[yb,yt] を [-0.9,0.9]×[-0.9,0.9] */
  xc = (xl + xr) / 2;
  yc = (yb + yt) / 2;
  ratio_x = 1.8 / (xr - xl);
  ratio_y = 1.8 / (yt - yb);
}

double X(double x) { return ratio_x*(x-xc); }
double Y(double y) { return ratio_y*(y-yc); }

void keyboard(unsigned char key, int x, int y)
{
  switch (key) {
  case 'q':
  case 'Q':
  case '\033':
    exit(0);
  default:
    break;
  }
}

void init(void)
{
  glClearColor(0.0, 0.0, 1.0, 1.0);
}

void resize(int w, int h)
{
  /* ウィンドウ全体をヴューポートにする */
  glViewport(0, 0, w, h);
  /* 変換行列の初期化 */
  glLoadIdentity();
  /* */
  glOrtho(- w / WIDTH, w / WIDTH, -h / HEIGHT, h / HEIGHT, -1.0, 1.0);
}

void drawBitmapString(double x, double y, void *font, char *str)
{
  glRasterPos2d(X(x), Y(y));
  while (*str)
    glutBitmapCharacter(font, *str++);
}

int N;
double lambda, Tmax;

void display(void)
{
  int i, n, nMax;
  double tau, h;
  double *u, *newu;
  double f(double);
  int win;
  char message[100];

  /* h, τ を計算する */
  h = 1.0 / N;
  tau = lambda * h * h;
  printf("τ=%g\n", tau);

  /* ベクトル 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);

  screen(-0.2, -0.2, 1.2, 1.2);

  glClear(GL_COLOR_BUFFER_BIT);

  /* 座標軸を表示する */
  glColor3d(0.0, 1.0, 0.0);
  glBegin(GL_LINES);
  glVertex2d(X(-0.1), Y(0.0));
  glVertex2d(X( 1.1), Y(0.0));
  glVertex2d(X(0.0), Y(-0.1));
  glVertex2d(X(0.0), Y( 1.1));
  glEnd();

  /* t=0 の状態を表示する */
  glColor3d(1.0, 0.0, 0.0);
  glBegin(GL_LINE_LOOP);
  glVertex2d(X(0.0), Y(u[0]));
  for (i = 1; i <= N; i++)
    glVertex2d(X(i * h), Y(u[i]));
  glEnd();

  /* タイトルと入力パラメーターを表示する */
#ifdef NONE
  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);
#endif
  drawBitmapString(0.1, 0.8, GLUT_BITMAP_TIMES_ROMAN_24,
		   "heat equation, u(0,t)=u(1,t)=0");

  /* ループを何回まわるか計算する (四捨五入) */
  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)τ) の状態を表示する */

    glBegin(GL_LINE_LOOP);    
    glVertex2d(X(0.0), Y(u[0]));
    for (i = 1; i <= N; i++)
      glVertex2d(X(i * h), Y(u[i]));
    glEnd();
  }
  glFlush();
}

int main(int argc, char **argv)
{
  /* N, λ を入力する */
  printf("区間の分割数 N = "); scanf("%d", &N);
  printf("λ (=τ/h^2) = "); scanf("%lf", &lambda);
  /* 最終時刻を入力する */
  printf("最終時刻 Tmax = "); scanf("%lf", &Tmax);
  printf("計算終了後、プログラムを終了するには q または ESC\n");
  /* ***************** グラフィックスの準備 ***************** */
  glutInitWindowPosition(0, 0);
  glutInitWindowSize(WIDTH, HEIGHT);
  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_RGBA);
  glutCreateWindow(argv[0]);
  glutDisplayFunc(display);
  glutReshapeFunc(resize);
  glutKeyboardFunc(keyboard);
  init();
  glutMainLoop();

}

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