3.7 プログラム

問題($A$)を解くためのプログラムである。
（これから載せるプログラムを実行するには、$N=500$以上取ることが お勧めである。それ以下だと波動はギザギザになる。）

/*
 * move1d-e.c -- 1次元波動方程式の初期値境界値問題を解く。
 *   コンパイル:  ccx move1d-e.c
 */

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

int main()
{
  int i, j, Jmax, skip, N;
  double tau, h, lambda, dt, Tmax;
  double *u, *newu, *v, *newv, *w, *neww;
  double f(double), g(double), fx(double); 

  /* N, λ を入力する */
  printf("区間の分割数 N = "); scanf("%d", &N);
  printf(" λ=  "); scanf("%lf", &lambda);
  
  /* hを計算する */
  h = 1.0 / N;
  tau = lambda * h;
  printf("τ=%f\n", tau);
  printf(" dt=  "); scanf("%lf", &dt);
  skip = dt / tau;

  /* 最終時刻を入力する */
  printf("最終時刻 Tmax = "); scanf("%lf", &Tmax);

  /* ベクトル u, newu, v, newv, w, neww を用意する */
  u = malloc(sizeof(double) * (N+1));
  newu = malloc(sizeof(double) * (N+1));
  v= malloc(sizeof(double) * (N+1));
  newv = malloc(sizeof(double) * (N+1));
  w= malloc(sizeof(double) * (N+1));
  neww = malloc(sizeof(double) * (N+1));

  for (i = 0; i <= N; i++)
    u[i] = f(i * h);

  for (i = 0; i <=N; i++)
    v[i] = g(i * h);

  for (i = 0; i <=N; i++)
    w[i] = fx(i * h);

  /* ウィンドウを出す */
  openpl();
  fspace2(-0.1, -1.1, 1.1, 1.1);

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

  /* ループを何回まわるか計算する */
  Jmax = Tmax / tau + 0.5;

  /* Dirichlet 境界条件 */
  u[0] = u[N] = 0.0;
  v[0] = v[N] = 0.0;

  for (j = 0; j <= Jmax; j++) {

    /* 差分方程式 (j -> j+1) */
    for (i = 1; i < N; i++){
      newv[i]= (v[i+1] + lambda* w[i+1])/ 2.0 
        +(v[i-1] - lambda* w[i-1])/ 2.0;
      neww[i]=( lambda * v[i+1] + w[i+1])/ 2.0
        +( -lambda * v[i-1] + w[i-1])/ 2.0;
      newu[i]=( u[i+1] + u[i-1])/ 2.0 + tau * v[i];
    } 

    neww[0] =  lambda * v[1] + w[1];
    neww[N] = -lambda * v[N-1] + w[N-1];

    for (i = 1; i < N; i++){
      u[i] = newu[i];
      v[i] = newv[i];
      w[i] = neww[i];
    }

    w[0] = neww[0];
    w[N] = neww[N];

    if (j % skip == 0) {
      /* この時刻 (t=(j+1)τ) の状態を表示する */
      erase();
      fmove(0.0, u[0]);
      for (i = 1; i <= N; i++)
        fcont(i * h, u[i]);
    }
  }
  printf("終りました。ウィンドウをクリックして下さい。\n");
  closepl();
}

double f(double x)
{
  if (x < 0.4)
    return 0.0;
  if (x > 0.6)
    return 0.0;
  else
    return sin(5.0 * M_PI * x);
}

double g(double x)
{
  return 0.0;
}

double  fx(double x)
{
  if (x < 0.4)
    return 0.0;
  if (x > 0.6)
    return 0.0;
  else
    return 5.0 * M_PI * cos(5.0 * M_PI * x);
}

問題($B$)を解くためのプログラムである。

/* move1d-e.c -- 1次元波動方程式の初期値境界値問題を解く。
 *   コンパイル:  ccx move1d-e.c
 * du/dx = 0 ( x= 0.0 )
 * du/dx = 0  ( x= 1.0 ) 
 *    の条件で解く。
 */

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

int main()
{
  int i, j, Jmax, skip, N;
  double tau, h, lambda, dt, Tmax, a, b, q, s;
  double *u, *newu, *v, *newv, *w, *neww;
  double f(double), g(double), fx(double); 

  /* N, λ を入力する */
  printf("区間の分割数 N = "); scanf("%d", &N);
  printf(" λ=  "); scanf("%lf", &lambda);

  /* hを計算する */
  h = 1.0 / N;
  tau = lambda * h;
  printf("τ=%f\n", tau);
  printf(" dt=  "); scanf("%lf", &dt);
  skip = dt / tau;

  /* 最終時刻を入力する */
  printf("最終時刻 Tmax = "); scanf("%lf", &Tmax);

  /* ベクトル u, newu, v, newv, w, neww を用意する */
  u = (double *)malloc(sizeof(double) * (N+1));
  newu = (double *)malloc(sizeof(double) * (N+1));
  v= (double *)malloc(sizeof(double) * (N+1));
  newv = (double *)malloc(sizeof(double) * (N+1));
  w= (double *)malloc(sizeof(double) * (N+1));
  neww = (double *)malloc(sizeof(double) * (N+1));

  for (i = 0; i <= N; i++)
    u[i] = f(i * h);

  for (i = 0; i <=N; i++)
    v[i] = g(i * h);

  for (i = 0; i <=N; i++)
    w[i] = fx(i * h);

  /* ウィンドウを出す */
  openpl();
  fspace2(-0.1, -1.1, 1.1, 1.1);

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

  /* ループを何回まわるか計算する */
  Jmax = (int)rint(Tmax / tau);

  /* Dirichlet 境界条件 */
  w[0] = 0.0;
  w[N] = 0.0;

  for (j = 0; j <= Jmax; j++) {

    /* 差分方程式 (j -> j+1) */
    for (i = 1; i < N; i++){
      newv[i]= (v[i+1] + lambda* w[i+1])/ 2.0 
        +(v[i-1] - lambda* w[i-1])/ 2.0;
      neww[i]=( lambda * v[i+1] + w[i+1])/ 2.0
        +( -lambda * v[i-1] + w[i-1])/ 2.0;
      newu[i]=( u[i+1] + u[i-1])/ 2.0 + tau * v[i];
    }

    newv[0]=  lambda * w[1] + v[1];
    newu[0]=  u[1] + tau * v[0];
    newv[N]= -lambda * w[N-1] + v[N-1];
    newu[N]=  u[N+1] + tau * v[N];
    neww[0]= w[0];
    neww[N]= w[N];

    for (i = 0; i <= N; i++){
      u[i] = newu[i];
      v[i] = newv[i];
      w[i] = neww[i];
    }

    if (j % skip == 0) {
    /* この時刻 (t=(j+1)τ) の状態を表示する */
      erase();
      fmove(0.0, u[0]);
      for (i = 1; i <= N; i++)
        fcont(i * h, u[i]);
    } 
  }
  printf("終りました。ウィンドウをクリックして下さい。\n");
  closepl();

}

double f(double x)
{
  if (x < 0.2)
    return 0.0;
  if (x > 0.6)
    return 0.0;
  else
    return sin(5.0 * M_PI * x);
 }

double g(double x)
{
  return 0.0;
}

double fx(double x)
{
  if (x < 0.2)
    return 0.0;
  if (x > 0.6)
    return 0.0;
  else
    return 5.0 * M_PI * cos(5.0 * M_PI * x);
}