Code:
//All Header Files Used
#include <iostream>
#include <cstdlib>
#include <string>
#include <cmath>
#include <fstream>
#include <cctype>
#include <climits>
//initialising all matrices and vectors to the size of 100
double A[100][100],L[100][100],f[100], d[100], X[100],U[100][100];
//declaring functions for use in the program
void Finite_Difference(int n,double a,double b,double c,double dx,double A[ ][100]);
void LU(int n, double A[][100],double L[][100]);
void f_sub(int n,double L[ ][100],double f[100],double g[100]);
void b_sub(int n,double U[ ][100],double g[100],double X[100]);
using namespace std;
//Function Definition to Calculate the Finite Difference Matrix
void Finite_Difference(int n,double a,double b,double c,double dx, double A[ ][100])
{
for(int i=0;i<n;i++)
{
for(int j=0;j<n;j++)
{
if (i == j)
{
A[i][j]=((2*a) - c*(pow(dx,2)));
}
else if( j == (i+1))
{
A[i][j] = (-a + (b/2)*dx);
}
else if (j == (i-1))
{
A[i][j] = (-a - (b/2)*dx);
}
else
A[i][j] = 0;
}
}
}
//Function Definition to Factor the A matrix in to the [L][U] matrices
void LU(int n,double U[ ][100],double L[ ][100] )
{
for(int j=0;j<(n-1);j++)
{
for(int i=j+1;i<n;i++)
{
L[i][j]= (U[i][j]/U[j][j]);
for(int k=j+1;k<n;k++)
{
U[i][k]=U[i][k]-(L[i][j]*U[j][k]);
}
}
}
for(int i=0;i<n;i++)
{
for(int j=0;j<n;j++)
{
if(i==j)
L[i][j]=1;
}
}
}
//Function Definition for Forward Substitution
void f_sub(int n,double L[ ][100],double f[100],double g[100] )
{
g[0]=f[0];
for(int i=1;i<n;i++)
{
g[i]=f[i];
for(int k=0;k<i;k++)
{
g[i]=g[i]-(L[i][k]*g[k] );
}
}
}
///Function Definition for Backward Substitution
void b_sub(int n,double U[ ][100],double g[100],double X[100] )
{
X[n-1]= g[n-1]/U[n-1][n-1];
for(int i=n-2;i>=0;i--)
{
X[i]=g[i];
for(int k=i+1;k<n;k++)
{
X[i]=X[i]-(U[i][k]*X[k]);
}
X[i]=X[i]/(U[i][i] );
}
}
int main()
{
//initialises a boolean so the user can choose to quit or continue another trial
bool again = true;
while(again)
{
cout<<" Welcome to the Boundary Value Problem Solver"<<endl;
cout<<" This Solver Will be Using the Finite Difference Method"<<endl;
cout<<" The Differential Equation to be Solved is in the Form of:"<<endl;
cout<<endl;
cout<<" ay'' + b' + cy + x = 0"<<endl;
cout<<" With Boundary Conditions y(0) = d & y(20) = e"<<endl;
double const1,const2,const3,bound1,bound2,dx;
bool v1 = false;
while(!v1)
{
cout<<"Please Enter the Constant 'a':\n";
cin>>const1;
//Determined whether the input is valid or invalid
if(cin.fail())
{
cin.clear();
cin.ignore(INT_MAX, '\n');
cout << "The input was invalid. Please try again" << endl;
}
else
{
v1 = true;
}
}
bool v2 = false;
while(!v2)
{
cout<<"Please Enter the Constant 'b':\n";
cin>>const2;
//Determined whether the input is valid or invalid
if(cin.fail())
{
cin.clear();
cin.ignore(INT_MAX, '\n');
cout << "The input was invalid. Please try again" << endl;
}
else
{
v2 = true;
}
}
bool v3 = false;
while(!v3)
{
cout<<"Please Enter the Constant 'c':\n";
cin>>const3;
//Determined whether the input is valid or invalid
if(cin.fail())
{
cin.clear();
cin.ignore(INT_MAX, '\n');
cout << "The input was invalid. Please try again" << endl;
}
else
{
v3 = true;
}
}
bool v4 = false;
while(!v4)
{
cout<<"Please Enter the Boundary Condtion at y(0):\n";
cin>>bound1;
//Determined whether the input is valid or invalid
if(cin.fail())
{
cin.clear();
cin.ignore(INT_MAX, '\n');
cout << "The input was invalid. Please try again" << endl;
}
else
{
v4 = true;
}
}
bool v5 = false;
while(!v5)
{
cout<<"Please Enter the Boundary Condition at y(20):\n";
cin>>bound2;
//Determined whether the input is valid or invalid
if(cin.fail())
{
cin.clear();
cin.ignore(INT_MAX, '\n');
cout << "The input was invalid. Please try again" << endl;
}
else
{
v5 = true;
}
}
bool v6 = false;
while(!v6)
{
cout<<"Please Enter the Interval Size:\n";
cin>>dx;
//Determined whether the input is valid or invalid
if(cin.fail())
{
cin.clear();
cin.ignore(INT_MAX, '\n');
cout << "The input was invalid. Please try again" << endl;
}
else if(dx<=0)
{
cout<< "The Interval Must be Greater than 0"<<endl;
cout<< "Please Try Again."<<endl;
}
else
{
v6 = true;
}
}
double dim = ((bound2-bound1)/dx) - 1;
int size = (int)dim;
//Initialises the A,L, and U matrices to 0 everywhere
for(int i=0;i<size-1;i++)
{
for(int j=0;j<size-1;j++)
{
A[i][j]=0;
L[i][j]=0;
U[i][j]=0;
}
}
//creates the Finite_Difference matrix
Finite_Difference(size,const1,const2,const3,dx,A);
//Create ouput text file
ofstream output_data;
output_data.open("output.txt");
//Check to see if output file has opened correctly
if(output_data.fail())
{
cout<<"Output file opening failed.\n";
exit(1);
}
//creates the Right-Hand Vector
int xi = 0;
for( i=0;i<size;i++)
{
xi+=dx;
if (i==0)
{
f[i] = bound1 + (xi*(pow(dx,2)));
}
else if (i == (size -1))
{
f[i] = bound2 +(xi*(pow(dx,2)));
}
else
{
f[i] = xi*(pow(dx,2));
}
}
cout<<endl;
LU(size,A, L);
f_sub( size,L, f, d);
b_sub( size,A, d, X);
output_data<<"The Solution is:"<<endl;
output_data<<" x y"<<endl;
output_data<<" _________ __________"<<endl;
output_data<<" "<<0<<" "<<bound1<<endl;
cout<<"The Solution to this Problem can be found in the file <output.txt>:\n";
int ex = dx;
for( i=0;i<size;i++)
{
ex+=dx;
output_data<<" "<<ex<<" "<<X[i]<<endl;
}
output_data<<" "<<20<<" "<<bound2<<endl;
cout<<"Would you like to Run This Program Again?(y/n)"<<endl;
string encore;
cin>>encore;
//restarts another analysis
if(encore == "y" || encore == "Y")
{
//clears the screen
system("cls");
again = true;
}
//quits the program
else if(encore == "n" || encore == "N")
{
again = false;
}
//Check whether the input is valid
else
{
cout<<"You Have not Entered a Proper Argument. Quiting Anyways"<<endl;
again = false;
}
}
return 0;
}
