Numeryczne rozwiązywanie układu równań liniowych metodą Jacobiego.

JacobianAndGaussSeidel.cpp

#include <iostream>
#include <cstdlib>
#include <time.h>
using namespace std;

//matrix of constatns:
//|a0, a1, a2, a3|
//|a4, a5, a6, a7|
//|a8, a9,a10,a11|
//matrix of variables: |x0,x1,x2...xn|
//x1 = (-a0x1 -a1x2-a3x3 +a0x1+a3)/a0
//x2 = (-a4x1 -a5x2-a6x3 +a5x2+a7)/a5
//x3 = (-a8x1 -a9x2-a10x3 +a10x3+a11)/a11

long double calculateVariableInRowOfMatrix( long double ** constantsMatrix, long double * variableValues, int numeberOfVariables, int indexOfVariable){
	long double valueOfXi=0;
	for (int i = 0; i < numeberOfVariables; i++){
		valueOfXi -= constantsMatrix[indexOfVariable][i]*variableValues[i];
	}
	valueOfXi += constantsMatrix[indexOfVariable][indexOfVariable] * variableValues[indexOfVariable];
	valueOfXi += constantsMatrix[indexOfVariable][numeberOfVariables];
	//cout << valueOfXi <<"divideby:" << constantsMatrix[indexOfVariable][indexOfVariable]<<endl;
	valueOfXi /= constantsMatrix[indexOfVariable][indexOfVariable];

	return valueOfXi;
}
void deepCopyItems1D(long double * &from, long double* &to,int size){
	for (int i = 0; i < size; i++){
		to[i] = from[i];
	}
}
long double * jacobianMethod(long double ** constantsMatrix,long double * initialVariableValues, int numberOfVaiables,int arbitraryNumberOfIterations){
	long double * previousItarationVariableValues = new long double[numberOfVaiables];
	long double * resultTable = new long double[numberOfVaiables];
	deepCopyItems1D(initialVariableValues, previousItarationVariableValues,numberOfVaiables);
	deepCopyItems1D(initialVariableValues, resultTable, numberOfVaiables);
	for (int i = 0; i < arbitraryNumberOfIterations; i++){
		for (int j = 0; j < numberOfVaiables; j++){
			resultTable[j] = calculateVariableInRowOfMatrix(constantsMatrix, previousItarationVariableValues, numberOfVaiables, j);
		}
		deepCopyItems1D(resultTable, previousItarationVariableValues, numberOfVaiables);
	}
	return resultTable;
}


long double * gaussSeidelMethod(long double ** constantsMatrix, long double * initialVariableValues, int numberOfVaiables, int arbitraryNumberOfIterations){
	long double * resultTable = new long double[numberOfVaiables];
	deepCopyItems1D(initialVariableValues, resultTable, numberOfVaiables);
	for (int i = 0; i < arbitraryNumberOfIterations; i++){
		for (int j = 0; j < numberOfVaiables; j++){
			resultTable[j] = calculateVariableInRowOfMatrix(constantsMatrix, resultTable, numberOfVaiables, j);
		}
	}
	return resultTable;
}



long double ** create2Dtable(int size){
	long double ** matrix = new long double *[size];
	for (int i = 0; i < size; i++)matrix[i] = new long double[size + 1];
	return matrix;
}
long double ** fillRandom2D(long double ** & matrix,int size){
	for (int i = 0; i < size; i++){
		for (int j = 0; j < size+1; j++){
			int n = (rand() % 51) - 25;
			if (n!=0)matrix[i][j] =n ;
			else j--;
		}
	}
	return matrix;
}
void print2D(long double ** matrix,int size){
	for (int i = 0; i < size; i++){
		for (int j = 0; j < size + 1; j++){
			cout << matrix[i][j]<<" ";
		}
		cout << endl;
	}
}
void print1D(long double * matrix, int size){
	for (int i = 0; i < size; i++){
			cout << matrix[i] << " ";
		}
}

void testJacobians(){
	int numberOfVariables = 3;
	long double ** matrix1 = create2Dtable(numberOfVariables);
	matrix1[0][0] = 5;			//x1
	matrix1[0][1] = -1;			//x2
	matrix1[0][2] = 2;			//x3
	matrix1[0][3] = 12;			//result

	matrix1[1][0] = 3;			//x1
	matrix1[1][1] = 8;			//x2
	matrix1[1][2] = -2;			//x3
	matrix1[1][3] = 25;			//result

	matrix1[2][0] = 1;			//x1
	matrix1[2][1] = 1;			//x2
	matrix1[2][2] = 4;			//x3
	matrix1[2][3] = 6;			//result
	//expected result: 
	//x1 = 2.7241379310344827586206896551724
	//x2 = 2.1724137931034482758620689655172
	//x3 = 0.27586206896551724137931034482759
	cout << "Given matrix:"<<endl;
	print2D(matrix1, numberOfVariables);
	long double initialVariableValues[3] = { 0, 0, 0 };
	cout << endl << "Jacobian method results: " << endl;
	print1D(jacobianMethod(matrix1, initialVariableValues, numberOfVariables, 100),numberOfVariables);
	cout << endl<< "Gauss Seidel method results: " << endl;
	print1D(gaussSeidelMethod(matrix1, initialVariableValues, numberOfVariables, 100), numberOfVariables);
	cout << endl;
}

void randomJacobians(int numberOfVariables, long double * initialVariableValues, int arbitraryNumberOfIterations){
	
	long double ** matrix1 = create2Dtable(numberOfVariables);
	matrix1 = fillRandom2D(matrix1, numberOfVariables);
			
	cout << endl << "Random matrix:" << endl;
	print2D(matrix1, numberOfVariables);
	cout << endl << "Jacobian method results: " << endl;
	print1D(jacobianMethod(matrix1, initialVariableValues, numberOfVariables, arbitraryNumberOfIterations), numberOfVariables);
	cout << endl << endl << "Gauss Seidel method results: " << endl;
	print1D(gaussSeidelMethod(matrix1, initialVariableValues, numberOfVariables, arbitraryNumberOfIterations), numberOfVariables);
	cout << endl;

}


int main(){
	srand(time(NULL));
	int numberOfVariables = 3;
	long double initialVariableValues[3] = { 0, 0, 0 };
	int arbitraryNumberOfIterations = 200; //works for <1000 because number becomes too long for long double (sometimes prints #IND when number has too much decimals)
	
	
	testJacobians();
	randomJacobians(numberOfVariables, initialVariableValues, arbitraryNumberOfIterations);
	randomJacobians(numberOfVariables, initialVariableValues, arbitraryNumberOfIterations);
	randomJacobians(numberOfVariables, initialVariableValues, arbitraryNumberOfIterations);
	randomJacobians(numberOfVariables, initialVariableValues, arbitraryNumberOfIterations);
	randomJacobians(numberOfVariables, initialVariableValues, arbitraryNumberOfIterations);

	system("pause");
}