///////////////////////////////////////////////////////////////////////////
//                                                                       //
// Program file name: Poisson.java                                       //
//                                                                       //
// © Tao Pang 2006                                                       //
//                                                                       //
// Last modified: January 18, 2006                                       //
//                                                                       //
// (1) This Java program is part of 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.     //
//                                                                       //
///////////////////////////////////////////////////////////////////////////

// Program for the one-dimensional Poisson equation.

import java.lang.*;
public class Poisson {
  final static int n = 99, m = 2;
  public static void main(String argv[]) {
  double d[] = new double[n];
  double b[] = new double[n];
  double c[] = new double[n];
  double h = 1.0/(n+1), pi = Math.PI;
  
 // Evaluate the coefficient matrix elements
    for (int i=0; i<n; ++i) {
      d[i] =  2;
      c[i] = -1;
      double xm = h*i;
      double xi = xm + h;
      double xp = xi + h;
      b[i] = 2*Math.sin(pi*xi) - Math.sin(pi*xm)
            -Math.sin(pi*xp);
    }

 // Obtain the solution
    double u[] = tridiagonalLinearEq(d, c, c, b);

 // Output the result
    double x = h;
    double mh = m*h;
    for (int i=0; i<n; i+=m) {
      System.out.println(x + " " + u[i]);
      x += mh;
    }
  }

// Method to solve the tridiagonal linear equation set.

  public static double[] tridiagonalLinearEq(double d[],
    double e[], double c[], double b[]) {
    int m = b.length;
    double w[] = new double[m];
    double y[] = new double[m];
    double z[] = new double[m];
    double v[] = new double[m-1];
    double u[] = new double[m-1];

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

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

 // Backward substitution to obtain z
    z[m-1] = y[m-1];
    for (int i=m-2; i>=0; --i) {
      z[i] = y[i]-u[i]*z[i+1];
    }
    return z;
  }
}