#include <stdio.h>
#include <stdlib.h>
#include <math.h>
void main()
{
int i,j,m,n
double TM=0,TMm=0,TN=0,TNn=0
int NN = 4
double L[4][4]= {{1,0,0,0}, {0,1,0,0}, {0,0,1,0},{0,0,0,1}, }
// M coef 4x4, N const 4x1
double M[4][4]={4,-1,0,2,-1,4,-1,0,0,-1,4,-1,2,0,-1,4}
double N[4][1]={-1,-7,9,0}
for(i=1i<NNi++) {
for(j=0j<ij++) {
for(m=0m<jm++) {
TM=M[i][m] * L[j][m] + TMm
TMm=TM
TM=0
}
M[i][j]=M[i][j]-TMm
TMm=0
L[i][j]=M[i][j]/M[j][j]
}
for(n=0n<in++) {
TN=M[i][n]*L[i][n]+TNn
TNn=TN
TN=0
}
M[i][i]=M[i][i]-TNn
TNn=0
}
// solve
for(i=0i<NNi++) {
for(m=0m<im++) {
TM=L[i][m]*N[m][0]+TMm
TMm=TM
TM=0
}
N[i][0]=N[i][0]-TMm
TMm=0
}
for(i=0i<NNi++) {
N[i][0]=N[i][0]/M[i][i]
}
for(i=NN-1i>=0i--) {
for(m=i+1m<NNm++) {
TM=L[m][i]*N[m][0]+TMm
TMm=TM
TM=0
}
N[i][0]=N[i][0]-TMm
TMm=0
}
for(i=0i<NNi++) printf("%0.2f\n",N[i][0])
}
这里向你推荐一下克鲁特算法（其实就是对高斯列主元消元法进行优化，使之更适合于计算机编程），首先将矩阵A进行LU分解（将系数矩阵分解成一个上三角矩阵和一个下三角矩阵），分解的过程中用到了隐式的主元寻找法，同时利用克鲁特算法可以将两个n*n矩阵压缩到一个n*n矩阵中，大大节省了存储空间提高了计算速度。方程可化为L*U*x=B，令U*x=y --->L*y=B
然后利用回代先求y，再利用y求x
因为该方法在求解过程中不涉及增广矩阵所以矩阵B几乎不参与什么运算，所以它的计算速度应该能够达到高斯列主元消元法的三倍，但原理与其基本一致。
而且我在程序中使用了动态数组方便你今后进行扩展。
以下程序按照《矩阵论第二版》和《C语言数值计算法方法大全》编写，LU分解部分程序主要参考了《C语言数值计算法方法大全》第二章的程序
如果你需要详细的理论讲解我可以将这两本书和源程序发给你，上面的论述相当详细足够你答辩用的了，我的邮箱[email protected]
计算结果：
A矩阵:
2 2 5
3 4 7
1 3 3
B矩阵:
5
6
5
解矩阵:
x 1=-7
x 2=0.333333
x 3=3.66667
Press any key to continue
#include <cmath>
#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <functional>
#include <vector>
#include <algorithm>
using namespace std
#define TINY 1.0e-20 //A small number.
#define N 3
void ludcmp(vector<vector<float>>&a, int n, vector<int>&indx, float &d)//对矩阵进行LU分解
void lubksb(vector<vector<float>>&a, int n, vector<int>&indx, vector<float>&b)//对矩阵进行前向和后向回代
void root(vector<vector<float>>&x,vector<float>&col)//解方程结果保存在y中
void iniv(vector<vector<float>>&x,vector<float>line,int n)//对二维动态数组进行初始化
void main()
{
int i,j,n=N//输入矩阵的维数
float A[N][N]={{2,2,5},{3,4,7},{1,3,3}}//左边A矩阵
float B[N]={5,6,5}//右边B矩阵
vector<vector<float>>x//建立动态二维数组存放A，保证你的程序进行扩展时只改A，B，N
vector<float>line
vector<float>y(n)//建立动态数组存放B
iniv(x,line,n)
y.clear()
for(i=0i<ni++)//将A赋给x，B赋给y
{
y.push_back(B[i])
for(j=0j<nj++)
{
x[i].push_back(A[i][j])
}
}
cout<<"A矩阵:"<<endl
for(i=0i<ni++)
{
for(j=0j<nj++)
{
cout<<setw(2)<<setiosflags(ios::left)<<setw(2)<<x[i][j]<<" "
}
cout<<endl
}
cout<<"B矩阵:"<<endl
for(i=0i<ni++)
{
cout<<setw(2)<<setiosflags(ios::left)<<setw(2)<<y[i]<<endl
}
root(x,y)//求根
cout<<"解矩阵:"<<endl
for(i=0i<ni++)
{
cout<<setw(2)<<setiosflags(ios::left)<<"x"<<i+1<<"="<<y[i]<<endl
}
cout<<endl
}
void root(vector<vector<float>>&x,vector<float>&col)
{
int n=x.size(),i=0,j=0
vector<int>index(n)//用于记录寻找主元素过程中对矩阵的初等变换
index.clear()
float m=1.0//记录变换方式，此程序中无用
ludcmp(x,n,index,m)//进行LU分解
lubksb(x,n,index,col)//根据分解结果进行回带
}
//以下程序按照《矩阵论第二版》和《C语言数值计算法方法大全》编写，LU分解部分程序主要参考了《C语言数值计算法方法大全》第二章的程序
//如果你需要详细的理论讲解我可以将这两本书和源程序发给你，我的邮箱[email protected]
void ludcmp(vector<vector<float>>&a, int n, vector<int>&indx, float &d)
{
int i,imax,j,k
float big=0,dum=0,sum=0,temp=0
vector<float>vv(n)
vv.clear()
d=1.0
for (i=0i<ni++)
{
big=0.0
for (j=0j<nj++)
if ((temp=fabs(a[i][j])) >big)
big=temp
vv[i]=1.0/big
}
for (j=0j<nj++)
{
for (i=0i<ji++)
{
sum=a[i][j]
for (k=0k<ik++)
sum -= a[i][k]*a[k][j]
a[i][j]=sum
}
big=0.0
for (i=ji<ni++)
{
sum=a[i][j]
for (k=0k<jk++)
sum -= a[i][k]*a[k][j]
a[i][j]=sum
if ( (dum=vv[i]*fabs(sum)) >= big)
{
big=dum
imax=i
}
}
if (j != imax)
{
for (k=0k<nk++)
{
dum=a[imax][k]
a[imax][k]=a[j][k]
a[j][k]=dum
}
d = -(d)
vv[imax]=vv[j]
}
indx[j]=imax
if (a[j][j] == 0.0)
a[j][j]=TINY
if (j != n)
{
dum=1.0/(a[j][j])
for (i=j+1i<ni++)
a[i][j] *= dum
}
}
}
void lubksb(vector<vector<float>>&a, int n, vector<int>&indx, vector<float>&b)
{
int i,ii=0,ip,j
float sum
for(i=0i<ni++)//按LU分解时寻找主元所进行的初等变换进行反边变换。
{
ip=indx[i]
sum=b[ip]
b[ip]=b[i]
b[i]=sum
}
sum=0
for (i=1i<ni++)
{
sum=0
for(j=0j<ij++)
{
sum+=a[i][j]*b[j]
}
b[i]=b[i]-sum
}
b[n-1]=b[n-1]/a[n-1][n-1]
for (i=n-2i>=0i--)
{
sum=0
for(j=i+1j<nj++)
{
sum+=a[i][j]*b[j]
}
b[i]=(b[i]-sum)/a[i][i]
}
}
void iniv(vector<vector<float>>&x,vector<float>line,int n)
{
int i,j
for(i=0i<ni++)
{
x.push_back(line)
for(j=0j<nj++)
{
x[i].clear()
}
}
}
