/**********
FullyImplicitDiffusion2: solves the diffusion equation for the test problem 
using the fully implicit difference method

This uses an entire array for the solution - compare Mathematica code uses a few vectors.

[in]: int M : number of time steps
      int N : number of positive and negative space steps

[out] : array of values in file for solution at given time

This uses Numerical Recipes Classes for Arrays
Numerical Recipes in C++, Second Edition

http://www.numerical-recipes.com/public_domain.html

Version 1.1nr, William Shaw, 2006.

*************/
#include <cmath>
#include <iostream>
#include <fstream>
#include "nrutil_nr.h" //allow access to public domain Numerical Recipes matrix routines

using namespace std;

void tridiagsolve(NRVec<double> &a, NRVec<double> &b, NRVec<double> &c, NRVec<double> &r, NRVec<double> &x)
{
	int j;
	int P=r.size();
	NRVec<double> u(P);
	NRVec<double> d(P);
	NRVec<double> l(P);
	NRVec<double> z(P);
	u[0]=c[0];
	d[0]=b[0];
	z[0]=r[0];
	// Do LU decomposition and forward subsititution as you go
	for (j=1;j<P;j++) {
		u[j]=c[j];
		l[j]=a[j]/d[j-1];
		d[j]=b[j]-l[j]*u[j-1];
		z[j] = r[j]-l[j]*z[j-1];
        }
	//follow with back substitution
	x[P-1] = z[P-1]/d[P-1];
	for (j=(P-2);j>=0;j--){
        x[j] = (z[j] - u[j]*x[j+1])/d[j];
        }
}



void ImplicitDiffusion(int N, int M)
{
       int i, j;
       double dt = 0.1/M ;
       double dx = 2.0/N;;
       const double PI = 3.141592653589793;
       NRMat<double> solvals(0.0,M+1,2*N+1); //NR dynamic matrix
       double alpha = 0.0; // key diffusion parameter
       NRVec<double> r(0.0,2*N-1); // create vectors for implicit solver
       NRVec<double> x(0.0,2*N-1);
       NRVec<double> a(0.0,2*N-1);
       NRVec<double> b(0.0,2*N-1);
       NRVec<double> c(0.0,2*N-1);

       alpha = dt/(dx*dx);
       
       //output dt
       cout << "dt = " << dt << "\n";
           //output dt
       cout << "dx = " << dx << "\n";
       // output alpha
       cout << "alpha = " << alpha << "\n";
       
       // initialize initital data
       
       for (j = 0; j <= 2*N; j++)
       {
           solvals[0][j] = sin(PI*(j-N)*dx/2);
       }
             
       // initialize boundary conditions
       
       for (i = 1; i< M; i++)
       {
           solvals[i][0] = 0;
           solvals[i][2*N] = 0;
           }
       // Run the fully implicit algorithm
       
       // First fill the arrays to feed to implicit solver
       // Note form for fully implicit scheme
       
       for (j=0; j<=2*N-2;j++)
       {
           a[j] = -alpha;
           b[j] = 1.0+2.0*alpha;
           c[j] = -alpha;
       }
       
       for (i= 1; i<= M; i++)
       {
           for (j = 1; j<= 2*N-1; j++)
           {
               r[j-1] = solvals[i-1][j];
               }// The boundary conditions are zero so no end effects
           tridiagsolve(a, b, c, r, x);
           
           for (j = 1; j<= 2*N-1; j++)
           {
           solvals[i][j] = x[j-1]; 
           }   
       }
       
    ofstream out("fdoutput3.txt");
    for (j=0; j<= 2*N; j++)
    {   out.precision(15);
        out << solvals[M][j] << "\n";
        }
       return ;
}
       
       
#include <iostream>

#include <cmath>
using namespace std;
int main()
{
    double result;
    char name[20];
    int N;
    int M;
    cout << "Implicit Diffusion Example \n ";
    cout << "Enter number of time steps (M) \n ";
    cin >> M;
    cout << M  << " time steps \n" ;
    cout << "Enter space step parameter (N) \n ";
    cin >> N;
    cout << 2*N  << " space steps \n" ;
    cout << " Result being written to fdoutput3.txt \n " << endl;
    ImplicitDiffusion(N,M);
    cout << "Hit any key+<RET> to finish \n ";
    cin >> name;
 
    return(0);
}