多元一次线性方程自然数解的算法
例如：
a1x1+a2x2+a3x3+a4x4=M
b1x1+b2x2+b3x3+b4x4=N
有满意加分
[解决办法]
这是我写的求逆矩阵的函数，剩下的就简单了。
#include <iostream>using namespace std;double* inverse_matrix(int determinant_scale, double *deter)//返回deter的逆矩阵,determinant_scale为矩阵的阶数，deter为第一个元素的地址{///////////////////////////////////////////////////////////////////////////////////////////初始化单位矩阵 double *i_matrix = new double[determinant_scale*determinant_scale]; for(int row_i=0; row_i<determinant_scale; ++row_i) { for(int column_i=0; column_i<determinant_scale; ++column_i) { if(row_i==column_i) { i_matrix[row_i*determinant_scale+column_i] = 1; } else { i_matrix[row_i*determinant_scale+column_i] = 0; } } }////////////////////////////////////////////////////////////////////////////////////////转变成上三角矩阵 for(int row=0; row<determinant_scale-1; ++row) { if(deter[row*determinant_scale+row]==0) { for(int row2=row+1; row2<determinant_scale && deter[row2*determinant_scale+row]==0; ++row2); double temp; for(int column2=row; column2<determinant_scale; ++column2) { temp = deter[row*determinant_scale+column2]; deter[row*determinant_scale+column2] = deter[row2*determinant_scale+column2]; deter[row2*determinant_scale+column2] = temp; } for(int column3=0; column3<determinant_scale; ++column3) { temp = i_matrix[row*determinant_scale+column3]; i_matrix[row*determinant_scale+column3] = i_matrix[row2*determinant_scale+column3]; i_matrix[row2*determinant_scale+column3] = temp; } } for(int row2=row+1; row2<determinant_scale; ++row2) { double k = -1*(deter[row2*determinant_scale+row]/deter[row*determinant_scale+row]); for(int column2=row; column2<determinant_scale; ++column2) { deter[row2*determinant_scale+column2] += k*deter[row*determinant_scale+column2]; } for(int column3=0; column3<determinant_scale; ++column3) { i_matrix[row2*determinant_scale+column3] += k*i_matrix[row*determinant_scale+column3]; } } }//////////////////////////////////////////////////////////////////////////////////////转变成对角矩阵 for(int re_row=determinant_scale-1; re_row>-1; --re_row) { for(int row2=re_row-1; row2>-1; --row2) { double k = -1*(deter[row2*determinant_scale+re_row]/deter[re_row*determinant_scale+re_row]); deter[row2*determinant_scale+re_row] += k*deter[row*determinant_scale+re_row]; for(int column3=0; column3<determinant_scale; ++column3) { i_matrix[row2*determinant_scale+column3] += k*i_matrix[re_row*determinant_scale+column3]; } } }/////////////////////////////////////////////////////////////////////////////////////////转变成单位矩阵 for(int r=0; r<determinant_scale; ++r) { for(int c=0; c<determinant_scale; ++c) { i_matrix[r*determinant_scale+c] *= (1/deter[r*determinant_scale+r]); } }///////////////////////////////////////////////////////////////////////////////////////////返回逆矩阵 return i_matrix;}int main(){ double a[3][3] = {1,2,3,2,1,2,1,3,3}; double *inverse_a; inverse_a = inverse_matrix(3, a[0]);//3为矩阵的阶数，a[0]为第一个元素的地址. for(int i=0; i<9; ++i) { cout<<inverse_a[i]<<" "; if(i%3==2) { cout<<endl; } } return 0;}
------解决方案--------------------
#include <iostream>#include <algorithm>#include <iterator>#include <vector>#include<windows.h>#include "..\LinearEquations\TMatrix.h"using namespace std;typedef double T;//从LU矩阵中获取L矩阵的元素T GetLElem(const TMatrix<T>& M, size_t x, size_t y){ if(x >= (M.GetXSize()) || y >= M.GetYSize()) throw range_error("下标错误！"); if (x == y) return 1; else if(x > y) return 0; else return M[y][x];}//从LU矩阵中获取U矩阵的元素T GetUElem(const TMatrix<T>& M, size_t x, size_t y){ if(x >= M.GetXSize() || y >= M.GetYSize()) throw range_error("下标错误！"); if(y > x ) return 0; else return M[y][x];}void MakeLU(TMatrix<T>& A){ //主对角线 for (size_t i = 0 ; i < (A.GetXSize() - 1) ; ++i) { //产生U矩阵的第i行 size_t& y = i; for (size_t x = y; x < A.GetXSize(); ++x) { T sum = 0; size_t t = 0; for (size_t xy = 0; xy < A.GetYSize(); ++xy) { if(xy == y) {//如果现在计算的是当前要生成的元素 t = xy ; continue; } sum += GetLElem(A, xy, y) * GetUElem(A, x, xy); } A[y][x] = (A[y][x] - sum)/* / GetLElem(A, t, y)*/; } //产生L矩阵的第i列 size_t& x = i; for (size_t y = x + 1; y < A.GetYSize(); ++y) { T sum = 0; size_t t = 0; for (size_t xy = 0 ; xy < A.GetYSize(); ++xy) { if(xy == x) { t = xy; continue; } sum += GetLElem(A, xy, y) * GetUElem(A, x, xy); } A[y][x] = (A[y][x] - sum) / GetUElem(A, x, t); } }}void SolveUpTrangle(TMatrix<T>& matrix, vector<T>& xs){ //从系数矩阵右下角的元素开始 for (int i = matrix.GetYSize() - 1; i >= 0 ; --i) { // T tmp = 0; for (int j = i + 1; j < matrix.GetXSize() - 1; j++) { tmp += matrix[i][j]/* * maxtrix[j][j]*/; } matrix[i][i] = (matrix[i].back() - tmp) / matrix[i][i]; for (int k = i - 1; k >= 0 ; --k) { matrix[k][i] *= matrix[i][i]; } } for (size_t i = 0 ; i < matrix.GetYSize() ; ++i) { xs.push_back(matrix[i][i]); }}void main(){ size_t N = 0; while (cout<<"请输入方程个数：", !(cin>>N)) { cout<<"数据错误，请重新输入！"<<endl; cin.clear(); cin.sync(); } TMatrix<T> A(N + 1, N); //TMatrix<T> LU(N + 1, N); A.In(); bool InputDone = false; do { cout<<"此方程组的增广矩阵为："<<endl; A.Out(); cout<<"请检查增广矩阵是否正确，正确则按回车继续，m键修改元素！"; char cmd; cin.clear(); cin.sync(); cin.get(cmd); switch(cmd) { case '\n':case '\r':case '\0': InputDone = true; break; case 'm':case 'M': unsigned x, y; while(cout<<"请先输入要修改元素的横、纵座标，Ctrl+z结束修改：",cin>>x>>y) { if(--x >= A.GetXSize() || --y >= A.GetYSize()) { cout<<"座标值错误，请重新输入!"<<endl; continue; } cin.sync(); cout<<"此位置原元素为："<<A[y][x] <<"，请输入新值："; while(!(cin>>A[y][x])) { cout<<"数据错误，请重新输入："; cin.clear(); cin.sync(); } } } } while (!InputDone); vector<T> Xs;// LARGE_INTEGER li1,li2,li;// memset(&li1, 0, sizeof(li));// memset(&li2, 0, sizeof(li));// memset(&li, 0, sizeof(li));// QueryPerformanceFrequency(&li);// QueryPerformanceCounter(&li1); MakeLU(A); SolveUpTrangle(A, Xs);// QueryPerformanceCounter(&li2); copy(Xs.begin(), Xs.end(), ostream_iterator<T>(cout, "\n"));// cout<<"求解用时："<<(li2.QuadPart - li1.QuadPart)/(double)li.QuadPart<<"秒"<<endl; system("pause");}