一，要解决的问题

选用合适的算法，求解三种线性方程组：一般线性方程组，对称正定方程组，三对角线性方程组。
方程略。

二，数值方法

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"<<" ";elsecout<<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];}//计算ydouble s3=0;for(int m=1;m<=k-1;m++)s3+=a[k][m]*yy[m];yy[k]=(b[k]-s3)/a[k][k];}//计算xxx[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"<<" ";elsecout<<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"<<" ";//elsecout<<res[i]<<" ";

能够看到x3的值是一个十分接近于0的数，假设将消元后的系数矩阵打印出来。能够看到消元后的系数矩阵并非一个真正的上三角矩阵，下三角部分有几处的值是一个绝对值极小的值。这是因为计算机的浮点计算造成的，浮点数在计算机中本身就不是一个精确的数，在消元的过程中。一些浮点运算有误差，于是最后得到的是近似值，而不是0。

同理，平方根法和追赶法也会产生由浮点数计算引起的误差，降低计算误差正是学习数值分析的目的。