一，要解决的问题
选用合适的算法，求解三种线性方程组：一般线性方程组，对称正定方程组，三对角线性方程组。
方程略。
二，数值方法
1，使用Guass列主元消去法求解一般线性方程组。
Guass列主元是为了防止Guass消去法中大数吃掉小数而引出的一种线性方程组求解方法，消元时选用一列中绝对值最大的元素作为列主元素。
算法伪代码：
消元过程
回代过程
2，使用平方根法求解对称正定方程组
平方根法。它把系数矩阵（对称正定矩阵）表示成一个下三角矩阵L和其转置的乘积的分解。这样的分解又称为Cholesky分解。
3，使用追赶法求解三对角线性方程组
三对角线性方程组是指这一类的线性方程组：系数矩阵是一个对角占优的三对角矩阵。追赶法是专门用来求解三对角线性方程组的，它将系数矩阵分解成alpha矩阵和beta矩阵的乘积，例如以下图所看到的：
三。算法
1。Guass列主元消去法
/*
CreateOn:2016/03/20
Author:linxiaobai
Function:linear equation solution
solution1:列主元Guass消去法求解一般线性方程组
*/
#include "stdafx.h"
# include<iostream>
# include<algorithm>
# include<fstream>
# include<iomanip>
# include<cmath>
using namespace std;
double a[15][15];
const int N=10;
double res[N+1];
void printArry(double a[][15])//打印增广矩阵
{
for(int i=1;i<=10;i++)
{
for(int j=1;j<=11;j++)
{
cout<<setw(3)<<a[i][j]<<" ";
}
cout<<endl;
}
}
int _tmain(int argc, _TCHAR* argv[])
{
cout<<"【运用列主元Guass求解一般线性方程组】"<<endl;
//读入增广矩阵
ifstream in;
in.open("data.txt");
if(!in)
{
cerr<<"file open failed!"<<endl;
return 0;
}
double x;
int i=1,j=1;
while(!in.eof())
{
in>>a[i][j];
j++;
if(j>11)
{
i++;j=1;
}
}
cout<<"增广矩阵："<<endl<<"++++++++++++++++++++++++++"<<endl;
printArry(a);
cout<<"++++++++++++++++++++++++++";
for(int k=1;k<=N-1;k++)
{
double tmp=abs(a[k][k]);
int ind=k;
for(int j=k;j<=N;j++)//找绝对值最大的行
{
if(abs(a[j][k])>tmp)
{tmp=abs(a[j][k]);ind=j;}
}
//若a[ind][k]=0,停止计算
if(a[ind][k]==0){cout<<"no unique solution"<<endl;return 0;}
//绝对值最大的行交换到第k行
if(ind!=k)
for(int j=1;j<=N+1;j++)
swap(a[ind][j],a[k][j]);
//消元计算
double p;
for(int i=k+1;i<=N;i++)
{
p=a[i][k]/a[k][k];
for(int j=k;j<=N+1;j++)
a[i][j]-=p*a[k][j];
}
}
if(a[N][N]==0)
{cout<<"no unique solution"<<endl;return 0;}
//回代求解
res[N]=a[N][N+1]/a[N][N];
double s;
for(int i=N-1;i>=1;i--)
{
s=0;
for(int j=i+1;j<=N;j++)
s+=a[i][j]*res[j];
res[i]=(a[i][N+1]-s)/a[i][i];
}
//输出解向量为
cout<<endl<<endl<<"解向量为："<<endl;
for(int i=1;i<=N;i++)
if(abs(res[i])<10e-14)
cout<<"0"<<" ";
else
cout<<res[i]<<" ";
cout<<endl;
return 0;
}
2，使用平方根算法解对称正定方程组
/*
CreateOn:2016/03/20
Author:linxiaobai
Function:linear equation solution
solution1:使用平方根算法解对称正定方程组
*/
# include"stdafx.h"
# include<iostream>
# include<fstream>
# include<cmath>
# include<iomanip>
using namespace std;
double a[10][10];
const int N=8;
double b[N+1],xx[N+1],yy[N+1];
void printArry(double a[][10])//输出系数矩阵
{
for(int i=1;i<=8;i++)
{
for(int j=1;j<=8;j++)
{
cout<<setw(3)<<a[i][j]<<" ";
}
cout<<endl;
}
}int main()
{
cout<<"【运用平方根算法解对称正定方程组】"<<endl;
/*读入系数矩阵*/
ifstream in;
in.open("data2.txt");
if(!in)
{
cerr<<"file open failed!"<<endl;
return 0;
}
int i=1,j=1;
while(i<=N)
{
in>>a[i][j];
j++;
if(j>8)
{
i++;j=1;
}
}
/*读入b[]*/
for(int i=1;i<=N;i++)
in>>b[i];
cout<<"系数矩阵："<<endl<<"++++++++++++++++++++++++++"<<endl;
printArry(a);
cout<<endl<<endl<<"b："<<endl;
for(int i=1;i<=N;i++)
cout<<setw(3)<<b[i]<<" ";
cout<<endl<<"++++++++++++++++++++++++++";
//開始计算。A=GG',G存在A的下三角
for(int k=1;k<=N;k++)
{
double s1=0;
for(int m=1;m<=k-1;m++)
s1+=pow(a[k][m],2);
a[k][k]=sqrt(a[k][k]-s1);
for(int i=k+1;i<=N;i++)
{
double s2=0;
for(int m=1;m<=k-1;m++)
s2+=a[i][m]*a[k][m];
a[i][k]=(a[i][k]-s2)/a[k][k];
}
//计算y
double s3=0;
for(int m=1;m<=k-1;m++)
s3+=a[k][m]*yy[m];
yy[k]=(b[k]-s3)/a[k][k];
}
//计算x
xx[N]=yy[N]/a[N][N];
for(int k=N-1;k>=1;k--)
{
double s4=0;
for(int m=k+1;m<=N;m++)
s4+=a[m][k]*xx[m];
xx[k]=(yy[k]-s4)/a[k][k];
}
cout<<endl<<"解向量为:"<<endl;
for(int i=1;i<=N;i++)
cout<<setw(3)<<xx[i]<<" ";
cout<<endl;
return 0;
}3，使用追赶法解三对角线性方程组
/*
CreateOn:2016/03/20
Author:linxiaobai
Function:linear equation solution
solution1:使用追赶法解三对角线性方程组
*/
# include"stdafx.h"
# include<iostream>
# include<fstream>
# include<iomanip>
# include<cmath>
using namespace std;
const int N=10;
double a[N+1],c[N+1],d[N+1],xx[N+1],yy[N+1];
int main()
{
cout<<"【运用追赶法解三对角线性方程组】"<<endl;
//a[N]:主对角线上的元素
for(int i=1;i<=N;i++)
a[i]=4;
//c[N]:上辅对角线上的元素
for(int i=1;i<=N-1;i++)
c[i]=-1;
//d[N]:下辅对角线上的元素
for(int i=2;i<=N;i++)
d[i]=-1;
double b[]={0,7,5,-13,2,6,-12,14,-4,5,-5};
cout<<"++++++++++++++++++++++++++"<<endl;
cout<<"主对角线上的元素a:"<<endl;
for(int i=1;i<=N;i++)
cout<<setw(3)<<a[i]<<" ";
cout<<endl<<endl<<"上辅对角线上的元素c:"<<endl;
for(int i=1;i<=N-1;i++)
cout<<setw(3)<<c[i]<<" ";
cout<<endl<<endl<<"下辅对角线上的元素d:"<<endl;
for(int i=2;i<=N;i++)
cout<<setw(3)<<d[i]<<" ";
cout<<endl<<endl<<"b:"<<endl;
for(int i=1;i<=N;i++)
cout<<setw(3)<<b[i]<<" ";
cout<<endl<<"++++++++++++++++++++++++++"<<endl;
/*開始计算*/
double alpha[N+1],beta[N+1];
alpha[1]=a[1];
for(int i=1;i<=N-1;i++)
{
beta[i]=c[i]/alpha[i];
alpha[i+1]=a[i+1]-d[i+1]*beta[i];
}
yy[1]=b[1]/alpha[1];
for(int i=2;i<=N;i++)
{
yy[i]=(b[i]-d[i]*yy[i-1])/alpha[i];
}
xx[N]=yy[N];
for(int i=N-1;i>=1;i--)
xx[i]=yy[i]-beta[i]*xx[i+1];
//解向量为：
cout<<endl<<"解向量为："<<endl;
for(int i=1;i<=N;i++)
if(abs(xx[i])<10e-14)
cout<<"0"<<" ";
else
cout<<xx[i]<<" ";
cout<<endl;
return 0;
}
四，数值结果
1，Guass列主元消去法
2，使用平方根算法解对称正定方程组
3，使用追赶法解三对角线性方程组
五，结果分析与实验总结
浮点计算产生的误差
在Guass消元算法之前的代码中，我使用了近似的方法，将绝对值小于10的-14次方的值近似为0，如今去掉这个处理。来看一下结果：
for(int i=1;i<=N;i++)
//if(abs(res[i])<10e-14)
//cout<<"0"<<" ";
//else
cout<<res[i]<<" ";
能够看到x3的值是一个十分接近于0的数，假设将消元后的系数矩阵打印出来。能够看到消元后的系数矩阵并非一个真正的上三角矩阵，下三角部分有几处的值是一个绝对值极小的值。这是因为计算机的浮点计算造成的，浮点数在计算机中本身就不是一个精确的数，在消元的过程中。一些浮点运算有误差，于是最后得到的是近似值，而不是0。
同理，平方根法和追赶法也会产生由浮点数计算引起的误差，降低计算误差正是学习数值分析的目的。