//		-*- Mode: Java -*-
// Author(s):	J. Robert Buchanan and Zhoude Shao
// Address:	Department of Mathematics
//		Millersville University
//      	P.O. Box 1002
//		Millersville, PA 17551-0302
// Phone:	717-871-7305
// FAX:		717-871-7948
// Email:	Robert.Buchanan@millersville.edu, Zhoude.Shao@millersville.edu
//
// Date:	10/07/2016
//
// Purpose:	Implicit finite difference method for solving the
//		homogeneous wave equation with Dirichlet boundary conditions
//		This program accompanies Example 14.3.
//
// To compile:	javac implicitwave.java
// To execute:	java implicitwave
//
import java.io.*;

public class implicitwave
{

    public static double f(double x)
    {
	return(x*(1.0-x))	;
    }	// end of function finitcond( )

    public static double g(double x)
    {
	return(Math.sin(5.0*Math.PI*x))	;
    }	// end of function g( )

    // Crout factorization for solving a tridiagonal system of equations
    public static void crout(
			     int nrows	,	// dimension of system
			     double [][] a,	// tridiagonal linear system
			     double [] x	// solution of system
			     )
    {
	int i;	// counter(s)
	double l[][] = new double[nrows][];	// lower triangular matrix
	double u[][] = new double[nrows][];	// upper triangular matrix
	double z[] = new double[nrows];
	for ( i=0 ; i<nrows ; i++ ) {
	    l[i] = new double[2];
	    u[i] = new double[2];
	}	// end of for ( i ) loop
	l[0][0] = a[0][0];
	u[0][1] = a[0][1]/l[0][0];
	z[0] = a[0][3]/l[0][0];
	for ( i=1 ; i<nrows-1 ; i++ ) {
	    l[i][0] = a[i][0];
	    l[i][1] = a[i][1]-l[i][0]*u[i-1][1];
	    u[i][1] = a[i][2]/l[i][1];
	    z[i] = (a[i][3]-l[i][0]*z[i-1])/l[i][1];
	}	// end of for ( i ) loop
	l[nrows-1][0] = a[nrows-1][1];
	l[nrows-1][1] = a[nrows-1][2]-l[nrows-1][0]*u[nrows-2][1];
	z[nrows-1] = (a[nrows-1][3]-l[nrows-1][0]*z[nrows-2])/l[nrows-1][1];
	x[nrows-1] = z[nrows-1];
	for ( i=nrows-2 ; i>-1 ; i-- ) {
	    x[i] = z[i]-u[i][1]*x[i+1];
	}	// end of for ( i ) loop
	return	;
    }	// end of function crout( )
    
    public static void main(String args[])
    {
	int i, j;	// counter(s)
	int nrows = 10-1;
	double h;	// spatial discretization step size
	double k;	// temporal discretization step size
	double r;	// method parameter
	double a[][] = new double[nrows][4];	// augmented tridiagonal matrix
	double b[][] = new double[nrows][4];	// augmented tridiagonal matrix
	double soln[][] = new double[3][nrows];	// result for new time

	// Set step sizes
	h = 1.0/((double)(nrows+1))	;
	k = h/2.0	;
	r = k/h	;
	// Set initial conditions and output
	System.out.printf("%8.4f ", 0.0)	;
	for ( j=0 ; j<nrows ; j++ ) {
	    soln[0][j] = f(((double)(j+1))*h)	;
	    System.out.printf("%8.4f ", soln[0][j])	;
	}	// end of for ( j ) loop
	System.out.printf("%8.4f\n", 0.0)	;
	// Populate augmented tridiagonal matrix
	// First row of matrix
	a[0][0] = 2.0*(2.0+r*r);
	a[0][1] = -r*r;
	a[0][2] = 0.0;
	b[0][0] = 4.0*(2.0-r*r);
	b[0][1] = 2.0*r*r;
	b[0][2] = 0.0;
	// Rows 2 through nrows-1 of matrix
	for ( i=1 ; i<(nrows-1) ; i++ ) {
	    a[i][0] = -r*r;
	    a[i][1] = 2.0*(2.0+r*r);
	    a[i][2] = -r*r;
	    b[i][0] = 2.0*r*r;
	    b[i][1] = 4.0*(2.0-r*r);
	    b[i][2] = 2.0*r*r;
	}	// end of for ( i ) loop
	// Last row of matrix
	a[nrows-1][0] = 0.0;
	a[nrows-1][1] = -r*r;
	a[nrows-1][2] = 2.0*(2.0+r*r);
	b[nrows-1][0] = 0.0;
	b[nrows-1][1] = 2.0*r*r;
	b[nrows-1][2] = 4.0*(2.0-r*r);
	// Determine the solution at t = k
	i = 0	;
	a[i][3] = 0.5*b[i][0]*f(((double)(i+1))*h)	;
	a[i][3]+= k*a[i][0]*g(((double)(i+1))*h)	;
	a[i][3]+= 0.5*b[i][1]*f(((double)(i+2))*h)	;
	a[i][3]+= k*a[i][1]*g(((double)(i+2))*h)	;
	for ( i=1 ; i<(nrows-1) ; i++ ) {
	    a[i][3] = 0.5*b[i][0]*f(((double)i)*h)	;
	    a[i][3]+= k*a[i][0]*g(((double)i)*h)	;
	    a[i][3]+= 0.5*b[i][1]*f(((double)(i+1))*h)	;
	    a[i][3]+= k*a[i][1]*g(((double)(i+1))*h)	;
	    a[i][3]+= 0.5*b[i][2]*f(((double)(i+2))*h)	;
	    a[i][3]+= k*a[i][2]*g(((double)(i+2))*h)	;
	}	// end of for ( i ) loop
	i = nrows - 1	;
	a[i][3] = 0.5*b[i][1]*f(((double)i)*h)	;
	a[i][3]+= k*a[i][1]*g(((double)i)*h)	;
	a[i][3]+= 0.5*b[i][2]*f(((double)(i+1))*h)	;
	a[i][3]+= k*a[i][2]*g(((double)(i+1))*h)	;
	crout(nrows,a,soln[1]);
	// Output the solution for t = k
	System.out.printf("%8.4f ", 0.0)	;
	for ( j=0 ; j<nrows ; j++ ) System.out.printf("%8.4f ", soln[1][j]);
	System.out.printf("%8.4f\n", 0.0)	;
	// Iterate for t = 2k, ..., 10k
	for ( j=2 ; j<10+1 ; j++ ) {
	    System.out.printf("%8.4f ",0.0);
	    i = 0	;
	    a[i][3] = b[i][0]*soln[(j-1)%3][i]	;
	    a[i][3]-= a[i][0]*soln[(j-2)%3][i]	;
	    a[i][3]+= b[i][1]*soln[(j-1)%3][i+1];
	    a[i][3]-= a[i][1]*soln[(j-2)%3][i+1];
	    for ( i=1 ; i<(nrows-1) ; i++ ) {
		a[i][3] = b[i][0]*soln[(j-1)%3][i-1]	;
		a[i][3]-= a[i][0]*soln[(j-2)%3][i-1]	;
		a[i][3]+= b[i][1]*soln[(j-1)%3][i]	;
		a[i][3]-= a[i][1]*soln[(j-2)%3][i]	;
		a[i][3]+= b[i][2]*soln[(j-1)%3][i+1]	;
		a[i][3]-= a[i][2]*soln[(j-2)%3][i+1]	;
	    }	// end of for ( i ) loop
	    i = nrows - 1	;
	    a[i][3] = b[i][1]*soln[(j-1)%3][i-1];
	    a[i][3]-= a[i][1]*soln[(j-2)%3][i-1];
	    a[i][3]+= b[i][2]*soln[(j-1)%3][i]	;
	    a[i][3]-= a[i][2]*soln[(j-2)%3][i]	;
	    crout(nrows,a,soln[j%3]);
	    for ( i=0 ; i<nrows ; i++ ) System.out.printf("%8.4f ", soln[j%3][i])	;
	    System.out.printf("%8.4f\n",0.0)	;
	}	// end of for ( j ) loop
    }	// end of function main( )
}

//
//	EOF
//