/*************************    Program 8.1    ****************************/
/*                                                                      */
/************************************************************************/
/* Please Note:                                                         */
/*                                                                      */
/* (1) This computer program is written by Tao Pang in conjunction with */
/*     his book, "An Introduction to Computational Physics," published  */
/*     by Cambridge University Press in 1997.                           */
/*                                                                      */
/* (2) No warranties, express or implied, are made for this program.    */
/*                                                                      */
/************************************************************************/

#include <stdio.h>
#include <math.h>
#define NMAX 99

main()
/* This program solves the one-dimensional Poisson equation
   with the Galerkin method as described in the text.
   Copyright (c) Tao Pang 1997. */
{
int i,n;
double pi,xl,h,d,e,b0,b1,xim,xi,xip;
double b[NMAX],a[NMAX],y[NMAX],w[NMAX],u[NMAX];

n = NMAX;
pi = 4*atan(1);
xl = 1;
h = xl/(n+1);
d = 2/h;
e = -1/h;
b0 = pi/h;
b1 = 1/h;

/* Find the elements in L and U */

w[0] = d;
u[0] = e/d;
for (i = 1; i < n; ++i)
  {
  w[i] = d-e*u[i-1];
  u[i] = e/w[i];
  }

/* Assign the array B */

for (i = 0; i < n; ++i)
  {
  xim = h*i;
  xi  = h*(i+1);
  xip = h*(i+2);
  b[i] = b0*cos(pi*xi)*(xim+xip-2*xi)
        +b1*(2*sin(pi*xi)-sin(pi*xim)-sin(pi*xip));
  }

/* Find the solution */

y[0] = b[0]/w[0];
for (i = 1; i < n; ++i)
  {
  y[i] = (b[i]-e*y[i-1])/w[i];
  }
a[n-1] = y[n-1];
for (i = n-2; i >= 0; i = i-1)
  {
  a[i] = y[i]-u[i]*a[i+1];
  printf("%16.8lf %16.8lf\n", h*(i+1),a[i]); 
  }
}