/*************************************************************************/
/* Program file name: poisson.c                                          */
/*                                                                       */
/* © Tao Pang 2006                                                       */
/*                                                                       */
/* Last modified: July 7, 2006                                           */
/*                                                                       */
/* (1) This C program is created for the book, "An Introduction to       */
/*     Computational Physics, 2nd Edition," written by Tao Pang and      */
/*     published by Cambridge University Press on January 19, 2006.      */
/*                                                                       */
/* (2) No warranties, express or implied, are made for this program.     */
/*                                                                       */
/*************************************************************************/

#include <stdio.h>
#include <math.h>
#define n 99
#define m 2

// Program for the one-dimensional Poisson equation.

main() {
  int i;
  double pi, h, xm, xi, xp;
  double d[n],b[n],c[n],u[n];
  void tridiagonalLinearEq();

  pi = 4*atan(1);
  h = 1.0/(n+1);

 // Evaluate the coefficient matrix elements
  for (i=0; i<n; ++i) {
    d[i] =  2;
    c[i] = -1;
    xm = h*i;
    xi = xm + h;
    xp = xi + h;
    b[i] = 2*sin(pi*xi) - sin(pi*xm) - sin(pi*xp);
  }

// Obtain the solution
  tridiagonalLinearEq(d, c, c, b, u);

// Output the result
  double x = h;
  double mh = m*h;
  for (i=0; i<n; i+=m) {
    printf("%16.8lf %16.8lf\n", x, u[i]); 
    x += mh;
  }
}

// Functione to solve the tridiagonal linear equation set.

void tridiagonalLinearEq(d, e, c, b, z)
  double d[n], e[n], c[n], b[n], z[n];
{
  int i;
  double y[n], w[n], v[n-1], t[n-1];

// Evaluate the elements in the LU decomposition
  w[0] = d[0];
  v[0]  = c[0];
  t[0]  = e[0]/w[0];
  for (i=1; i<n-1; ++i) {
    w[i]  = d[i]-v[i-1]*t[i-1];
    v[i]  = c[i];
    t[i]  = e[i]/w[i];
  }
  w[n-1]  = d[n-1]-v[n-2]*t[n-2];

// Forward substitution to obtain y
  y[0] = b[0]/w[0];
  for (i=1; i<n; ++i)
    y[i] = (b[i]-v[i-1]*y[i-1])/w[i];

// Backward substitution to obtain z
  z[n-1] = y[n-1];
  for (i=n-2; i>=0; --i) z[i] = y[i]-t[i]*z[i+1];
}