// LinearEqu.cpp: implementation of the LinearEqu class.
//
//#include "stdafx.h"
#include "LinearEqu.h"
#include <iostream>
#include <cmath>
#include <vector>
using namespace std;//
// Construction/Destruction
//LinearEqu::LinearEqu(int size) : Matrix(size)
{sums = new double[size];solution = new double[size];
}LinearEqu::~LinearEqu()
{delete[] sums;delete[] solution;
}void LinearEqu::setLinearEqu(const double* a, const double* b){setMatrix(a);for(int i = 0; i < this->getSize(); i++)sums[i] = b[i];
}void LinearEqu::displayLinearEqu() const{cout << "The Linear equation is:" << endl;for(int i = 0; i < this->getSize(); i++){for(int j = 0; j < this->getSize(); j++)cout << element(i, j) << "\t";cout << sums[i] << endl;}
}void LinearEqu::displaySolution() const{cout << "The Result is:" << endl;for(int i = 0; i < this->getSize(); i++)cout << "x[" << i << "] = " << solution[i] << endl;
}//交换
inline void swap(double &v1, double &v2){double temp = v1;v1 = v2;v2 = temp;
}//列主元消去法
bool LinearEqu::solve(){int n = this->getSize();for(int k = 0; k < n - 1; k++){int is;double max = 0;for(int i = k, j = k; i < n; i++){double t = fabs(element(i, j));if(t > max){max = t;is = i;}}if(max == 0){return false;}//交换行else{if(is != k){for(int j = k; j < n; j++)swap(element(k, j), element(is, j));swap(sums[k], sums[is]);}}//消去double major = element(k, k);for(i = k + 1; i < n; i++){double m = element(i, k) / major;for(int j = k; j < n; j++){element(i, j) -= element(k, j) * m;}sums[i] -= sums[k] * m;}}//判断剩下的一个元素是否为0double d = element(n - 1, n - 1);if(fabs(d) < 1e-15){return false;}//回代过程solution[n - 1] = sums[n - 1] / d;for(int i = n - 2; i >= 0; i--){double t = 0.0;for(int j = i + 1; j <= n - 1; j++)t += element(i, j) * solution[j];solution[i] = (sums[i] - t)/element(i, i);}return true;
} //LU分解法
bool LinearEqu::LU(){int n = this->getSize();vector<vector <double> > L(n ,vector<double>(n));vector<vector <double> > U(n ,vector<double>(n));//U第一行赋值for(int i = 0; i < n; i++)U[0][i] = element(0, i);//L第一列赋值for(int j = 1; j < n; j++)L[j][0] = element(j, 0) / U[0][0];//剩余赋值for(i = 1; i < n; i++){for(j = i; j < n; j++){double sum1 = 0;for(int k = 0; k < i; k++)sum1 += L[i][k] * U[k][j];U[i][j] = element(i, j) - sum1;}for(j = i+1; j < n; j++){double sum2 = 0;for(int k = 0; k < i; k++)sum2 += L[j][k] * U[k][i];L[j][i] = (element(j, i) - sum2) / U[i][i];}}vector<double> Y(n);Y[0] = sums[0];for (i = 1; i<n; i++){double sum3 = 0;for (int k = 0; k<i; k++)sum3 += L[i][k] * Y[k];Y[i] = sums[i] - sum3;}solution[n - 1] = Y[n - 1] / U[n - 1][n - 1];for (i = n - 2; i >= 0; i--){double sum4 = 0;for (int k = i + 1; k<n; k++)sum4 += U[i][k] * solution[k];solution[i] = (Y[i] - sum4) / U[i][i];}return true;
}//雅可比迭代法
bool LinearEqu::Jacobi()
{int n = this->getSize();int i,j;vector<double> x1(n);vector<double> x2(n);vector<double> b(n);vector<vector <double> > a(n ,vector<double>(n));double sum, max=INT_MAX, eps=1e-5, count=0;for(i = 0; i < n; i++){b[i] = sums[i] / element(i, i);for(j = 0; j < n; j++){a[i][j] = - element(i, j) / element(i, i);}a[i][i] = 0;}while(max > eps){for(i = 0; i < n; i++){sum = 0;for(j = 0; j < n; j++)sum += a[i][j] * x1[j];x2[i]  = sum + b[i];solution[i] = x2[i];}max = fabs(x1[0] - x2[0]);for(i = 1; i < n; i++)if(max < fabs(x1[i] - x2[i]))max = fabs(x1[i] - x2[i]);for(i = 0; i < n; i++){x1[i] = x2[i];}count++;}cout << "经过" << count << "次迭代" << endl;return true;
}//高斯塞德尔迭代法
bool LinearEqu::GS()
{int n = this->getSize();int i,j,k;vector<double> x1(n);vector<double> x2(n);vector<double> b(n);vector<vector <double> > a(n ,vector<double>(n));double sum, max=INT_MAX, eps=1e-5, count=0;for(i = 0; i < n; i++){b[i] = sums[i] / element(i, i);for(j = 0; j < n; j++){a[i][j] = - element(i, j) / element(i, i);}a[i][i] = 0;}while(max > eps){for(i = 0; i < n; i++){sum = 0;for(j = 0; j < i; j++){sum += a[i][j] * x2[j];}for(k = j; k < n; k++){sum += a[i][k] * x1[k];}x2[i]  = sum + b[i];solution[i] = x2[i];}max = fabs(x1[0] - x2[0]);for(i = 1; i < n; i++)if(max < fabs(x1[i] - x2[i]))max = fabs(x1[i] - x2[i]);for(i = 0; i < n; i++){x1[i] = x2[i];}count++;}cout << "经过" << count << "次迭代" << endl;return true;
}