/*
 * heat1g.c 1次元熱方程式, 境界条件は Robin と Dirichlet 同時に扱う
 *
 */

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

double exact(double x, double t);
void trilu(int n, double *al, double *ad, double *au);
void trisol(int n, double *al, double *ad, double *au, double *b);

/* left, right は != 0 だと
 *
 *  ∂u
 *  ---- + h u=T
 *  ∂n
 *
 * を課す。== 0 だと
 *
 *  u=U
 *
 * を課す。
 */

int left = 0, right = 0;
double pi;

int main()
{
  int i,Nx,n,nMax;
  double *u, *uu;
  double *ad,*al,*au;
  double x,dx,lambda,dt,h,t,error,maxerror;
  double T0,T1,U0,U1;
  double Tmax = 0.1;

  pi = 4 * atan(1.0);
  h = 1;

  if (left)
    T0 = 0;
  else
    U0 = 0;
  if (right)
    T1 = 0;
  else
    U1 = 0;

  printf("Nx="); scanf("%d", &Nx);
  dx = 1.0 / Nx;
  lambda = 0.5;
  dt = lambda * dx * dx;
  nMax = Tmax / dt + 0.5;

  u = malloc(sizeof(double) * (Nx+1));
  uu = malloc(sizeof(double) * (Nx+1));
  ad = malloc(sizeof(double) * (Nx+1));
  au = malloc(sizeof(double) * (Nx+1));
  al = malloc(sizeof(double) * (Nx+1));

  /* 係数行列を作る */
  for (i = 1; i < Nx; i++) {
    ad[i] = 1+2*lambda;
    au[i] = al[i] = - lambda;
  }

  if (left) {
    /* 左端を周囲温度指定するため */
    ad[0] = 1 + 2*lambda*(1+h*dx);
    au[0] = - 2 * lambda;
  }
  else {
    ad[0] = 1; au[0] = 0;
  }
  if (right) {
    /* 右端を周囲温度指定するため */
    ad[Nx] = 1 + 2*lambda*(1+h*dx);
    al[Nx] = -2 * lambda;
  }
  else {
    ad[Nx] = 1; al[Nx] = 0;
  }
  trilu(Nx + 1, al, ad, au);

  /* 初期条件 */
  for (i = 0; i <= Nx; i++)
    uu[i] = u[i] = exact(i * dx, 0.0);

  for (n = 1; n <= nMax; n++) {
    /* uu には u と同じ値が入っている。
     * 最初と最後だけいじると連立1次方程式の右辺になる
     */
    if (left)
      uu[0] += 2 * dx * lambda * T0;
    else
      uu[0] = U0;
    if (right)
      uu[Nx] += 2 * dx * lambda * T1;
    else
      uu[Nx] = U1;

    trisol(Nx + 1, al, ad, au, uu);
    
    for (i = 0; i <= Nx; i++)
      u[i] = uu[i];
    t = n * dt;
#ifdef VERBOSE
    printf("%f\n", t);
    for (i = 0; i <= Nx; i++)
      printf("%5.2f ", u[i]);
    printf("\n");
    for (i = 0; i <= Nx; i++)
      printf("%5.2f ", exact(i * dx, t));
    printf("\n");
#endif
    maxerror = 0;
    for (i = 0; i <= Nx; i++) {
      error = fabs(exact(i * dx, t) - u[i]);
      if (error > maxerror)
	maxerror = error;
    }
    printf("error=%e\n", maxerror);
  }
}

#define lambda1 2.028757838110434
#define lambda2 4.913180439434883

double exact(double x, double t)
{
  if (left) {
    /* まだ厳密解を知りません */
  }
  else {
    if (right)
      return sin(lambda1 * x) * exp(- lambda1 * lambda1 * t)
	+ sin(lambda2 * x) * exp(- lambda2 * lambda2 * t);
    else
      return sin(pi * x) * exp(- pi * pi * t)
	+ sin(2 * pi * x) * exp(- 4 * pi * pi * t);
  }

}

/* 3項方程式 (係数行列が三重対角である連立1次方程式のこと) Ax=b を解く
 *
 *   入力
 *     n: 未知数の個数
 *     al,ad,au: 連立1次方程式の係数行列
 *       (al: 対角線の下側 i.e. 下三角部分 (lower part)
 *        ad: 対角線       i.e. 対角部分 (diagonal part)
 *        au: 対角線の上側 i.e. 上三角部分 (upper part)
 *        つまり
 *
 *         ad[0] au[0]   0   .................. 0
 *         al[1] ad[1] au[1]   0  ............. 0
 *           0   al[2] ad[2] au[2]  0 ......... 0
 *                         ....................
 *                         al[n-2] ad[n-2] au[n-2]
 *                           0     al[n-1] ad[n-1]

 *        al[i] = A_{i,i-1}, ad[i] = A_{i,i}, au[i] = A_{i,i+1},
 *        al[0], au[n-1] は意味がない)
 *
 *     b: 連立1次方程式の右辺の既知ベクトル
 *        (添字は 0 から。i.e. b[0],b[1],...,b[n-1] にデータが入っている。)
 *   出力
 *     al,ad,au: 入力した係数行列を LU 分解したもの
 *     b: 連立1次方程式の解
 *   能書き
 *     一度 call すると係数行列を LU 分解したものが返されるので、
 *     以後は同じ係数行列に関する連立1次方程式を解くために、
 *     関数 trisol() が使える。
 *   注意
 *     Gauss の消去法を用いているが、ピボットの選択等はしていな
 *     いので、ピボットの選択をしていないので、係数行列が正定値である
 *     などの適切な条件がない場合は結果が保証できない。
 */

/* 三重対角行列の LU 分解 (pivoting なし) */
void trilu(int n, double *al, double *ad, double *au)
{
  int i, nm1 = n - 1;
  /* 前進消去 (forward elimination) */
  for (i = 0; i < nm1; i++) {
    al[i + 1] /= ad[i];
    ad[i + 1] -= au[i] * al[i + 1];
  }
}

/* LU 分解済みの三重対角行列を係数に持つ3項方程式を解く */
void trisol(int n, double *al, double *ad, double *au, double *b)
{
  int i, nm1 = n - 1;
  /* 前進消去 (forward elimination) */
  for (i = 0; i < nm1; i++) b[i + 1] -= b[i] * al[i + 1];
  /* 後退代入 (backward substitution) */
  b[nm1] /= ad[nm1];
  for (i = n - 2; i >= 0; i--) b[i] = (b[i] - au[i] * b[i + 1]) / ad[i];
}