# Initializing and Solving matrices

• 11-16-2003
scottmanc
Initializing and Solving matrices
here is my whole program. It will be used to solve bvp's (boundary Value Problems) using the finite difference method. The problem i have is initializing the matrix that needs to be solved and then having it printed out. I've spent a lot of time on it today and i hope to figure it out later on. If anyone has any suggestions on how to improve the code so that i can actually get the matrix initialized and the solution printed out, all help will be greatly appreciated. By the way, the formulas used are correct
and the only thing that is having problems is initializing the A matrix and printing out the X matrix, which is the solution.

Thanks

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; }```
• 11-16-2003
scottmanc
oh and the method of solving the matrix used in the program is by LU decomposition.