Initializing and Solving matrices

This is a discussion on Initializing and Solving matrices within the C++ Programming forums, part of the General Programming Boards category; here is my whole program. It will be used to solve bvp's (boundary Value Problems) using the finite difference method. ...

  1. #1
    Registered User
    Join Date
    Feb 2003
    Posts
    60

    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;
    }
    C++ can hurt.

  2. #2
    Registered User
    Join Date
    Feb 2003
    Posts
    60
    oh and the method of solving the matrix used in the program is by LU decomposition.
    C++ can hurt.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. trouble initializing constructor
    By dantestwin in forum C++ Programming
    Replies: 11
    Last Post: 07-06-2004, 02:31 AM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21