/*
 * heat1d-e-ng.c -- 1次元熱伝導方程式の初期値境界値問題
 * http://nalab.mind.meiji.ac.jp/~mk/program/fdm/heat1d-e-ng.c
 *
 *  u_t(x,t)=u_{xx}(x,t)        x∈(0,1), t∈(0,∞)
 *  u(0,t)=u(1,t)=0             t∈(0,∞)
 *  u(x,0)=f(x)                 x∈[0,1]
 *
 * を陽解法で解く。グラフィックスもなしの、一番シンプルな C プログラム
 *
 *   コンパイル:  gcc heat1d-e-ng.c -lm
 *   実行:        ./a.out
 *   入力データ例: N=10, λ=0.5, Tmax=0.1 
 */

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

#define MAXN (100)

double u[MAXN+1], newu[MAXN+1];

void print_data(double t, int N, double u[])
{
  int i;
  printf("t=%f\n", t);
  for (i = 0; i <= N; i++)
    printf("%f ", u[i]);
  printf("\n");
}

int main()
{
  int i, n, nMax, N;
  double tau, h, lambda, Tmax;
  double f(double);

  /* 分割数 N, λ(=τ/h^2) を入力する */
  printf("区間の分割数 N (≦%d): ", MAXN); scanf("%d", &N);
  if (N > MAXN) {
    fprintf(stderr, "N が大きすぎます。\n");
    exit(1);
  }
  printf("λ (=τ/h^2): "); scanf("%lf", &lambda);

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

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

  /* 初期値の代入 */
  for (i = 0; i <= N; i++)
    u[i] = f(i * h);

  /* t=0 の状態を表示する */
  print_data(0.0, N, u);

  /* ループを何回まわるか計算する (四捨五入) */
  nMax = (int)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)τ) の状態を表示する */
    print_data((n + 1) * tau, N, u);
  }
  return 0;
}

/* 初期条件 --- ここでは min(x,1-x)=1/2-|x-1/2| という中央が高い山形の関数 */
double f(double x)
{
  if (x <= 0.5)
    return x;
  else
    return 1.0 - x;
}