1.4 wave1d.c

/*
 * wave1d.c --- 1次元波動方程式 (Dirichlet 境界条件)
 *
 *   ccx wave1d.c -I/usr/local/include -lmatrix
 * ノートパソコンでは
 *   ccx wave1d.c -I/usr/local/include -L/usr/local/lib -lmatrix
 */

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

void copy_vector(int, vector, vector);

double pi;

int main()
{
  /* N は区間の分割数 */
  int N;
  int i, n, nmax, nfunc;
  /* lambda はλ, tau はτ */
  double lambda, Tmax, h, tau, lambda2;
  /* u1[i]=u_{i,n-1}
   * u2[i]=u_{i,n}
   * u3[i]=u_{i,n+1]
   */
  vector u1, u2, u3;
  double phi(double, int), psi(double, int);

  pi = 4.0 * atan(1.0);

  /* N は 100 とか 1000 とか */
  printf("N="); scanf("%d", &N);
  u1 = new_vector(N+1);
  u2 = new_vector(N+1);
  u3 = new_vector(N+1);

  /* λについては本に詳しいが、λ=1 が良い結果を出す */
  printf("λ="); scanf("%lf", &lambda);
  /* Tmax は最終時刻だから、長く計算したければ大きな値を入れる */
  printf("Tmax="); scanf("%lf", &Tmax);

  lambda2 = lambda * lambda;
  h = 1.0 / N;
  tau = lambda * h;

  /* 初期値を選ぶ */
  printf("nfunc(初期値の種類 1..2)="); scanf("%d", &nfunc);
  /* 初期値 */
  /* u_{i,0}=φi */
  for (i = 0; i <= N; i++)
    u1[i] = phi(i * h, nfunc);
  /* u_{i,1}=(1-λ2)... */
  for (i = 1; i < N; i++)
    u2[i] = (1-lambda2) * u1[i]
      + 0.5 * lambda2 * (u1[i-1] + u1[i+1]) + tau * psi(i * h, nfunc);
  u2[0] = u2[N] = 0.0;

  /* グラフを描くための準備 */
  openpl();
  fspace2(-0.2, -1.2, 1.2, 1.2);

  /* t=0 でのグラフ */
  fmove(0.0, u1[0]);
  for (i = 1; i <= N; i++)
    fcont(i * h, u1[i]);

  /* t=t1=τ でのグラフ */
  erase();
  fmove(0.0, u1[0]);
  for (i = 1; i <= N; i++)
    fcont(i * h, u2[i]);

  nmax = (int)rint(Tmax / tau);
  /* 時間に関するループ */
  for (n = 1; n <= nmax; n++) {
    for (i = 1; i < N; i++)
      u3[i] = 2 * (1 - lambda2) * u2[i] + lambda2 * (u2[i-1] + u2[i+1]) - u1[i];
    u3[0] = u3[N] = 0.0;

    /* t=t_{n+1}=(n+1)τ でのグラフ */
    erase();
    fmove(0.0, u3[0]);
    for (i = 1; i <= N; i++)
      fcont(i * h, u3[i]);

    /* u1 <- u2, u2 <- u3 */
    copy_vector(N + 1, u1, u2);
    copy_vector(N + 1, u2, u3);
  }
  closepl();

  return 0.0;
}

/* φ */
double phi(double x, int nfunc)
{
  if (nfunc == 1) {
    return sin(pi * x);
  }
  else if (nfunc == 2) {
    if (x > 0.375 && x <= 0.5)
      return 8.0 * x - 3.0;
    else if (x > 0.5 && x < 0.625)
      return - 8.0 * x + 5.0;
    else
      return 0.0;
  }
  else
    return 0.0;
}

/* ψ */
double psi(double x, int nfunc)
{
  if (nfunc == 1) {
    return 0.0;
  }
  else if (nfunc == 2) {
    return 0.0;
  }
  else
    return 0.0;
}

/*
 * ベクトル b をベクトル a にコピー
 */
void copy_vector(int N, vector a, vector b)
{
  int i;
  for (i = 0; i < N; i++)
    a[i] = b[i];
  return;
}