病态线性方程的求解

科学与工程计算实验报告学号：姓名：1004111202 王巧1004111204 梅亚 1004111208 李兵一、实验内容：考虑方程组Hx=b 的求解，其中系数矩阵H 为Hilbert 矩阵，nj i j i h h H j i n n j i ,,2,1,,11,)(,, =-+==⨯这是一个著名的病态问题。

通过首先给定解（例如取为各个分量均为1）再计算出右端b 的办法给出确定的问题。

实验要求： （1）选择问题的维数为6，分别用Jacobi 迭代法、GS 迭代法和SOR 迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比较，结论如何？（2）逐步增大问题的维数，仍然用上述的方法来解它们，计算的结果如何？计算的结果说明了什么？（3）讨论病态问题求解的算法。

二、程序设计的基本思想、原理和算法描述：1、 算法 Jacobi 迭代法若A 为稀疏矩阵，只需遍历非零元素 GS 迭代法若A 为稀疏矩阵，只需遍历非零元素每步迭代计算量相当于一次矩阵与向量的乘法；不需要保留上一步的迭代解，与Jacobi 迭代法计算量一样。

SOR 迭代法（稠密矩阵）2、 函数组成double max(double array[100]) 求数组中的最大值函数 3、 输入/输出设计对于方程组Hx=b 的求解，系数矩阵H 为Hilbert 矩阵，矩阵中的数由下列函数生成。

nj i j i h h H j i n n j i ,,2,1,,11,)(,, =-+==⨯X*取一个特解[1,1,1, (1)b 数组由矩阵的每行元素相加得来。

4、符号名说明double c[100] 用来存储第k+1次和第k 次迭代结果的差值的绝对值 double x[100] 第k+1次迭代结果，储存解数组 double x0[100] 初始向量double r 第k+1次和第k 次迭代结果的差值的绝对值的最大值 double sum 矩阵方程变换后右侧值的和 int k 迭代次数double a[100][100] 存储Hilbert 矩阵 double b[100] 存储b 向量三、源程序及注释：Jacobi 迭代法 #include<iostream> #include<math.h> #include <iomanip> #include <stdio.h> using namespace std; int n;double max(double array[100])//求最大值函数 {double a=array[0]; int i;for(i=1; i<n; i++) {if(a<array[i]) a=array[i]; return a; } }double c[100]= {0.0};double x[100]= {0.0}; //第k+1次迭代结果，储存解数组 double x0[100]= {0.0}; //初始向量double r,sum=0;int main(){double s=0;int i,k,j;//k为迭代次数double a[100][100];double b[100]= {0.0};cout<<"请输入维数："<<endl;cin>>n;cout<<"输出a数组："<<endl;double max(double array[100]);for(i=0; i<n; i++){for(j=0; j<n; j++){a[i][j]=1.0/(i+j+1);printf("%10.6lf ",a[i][j]);//矩阵中的数精确到六位}cout<<endl;}cout<<"输出b数组："<<endl;for(i=0; i<n; i++){for(j=0; j<n; j++){b[i]+=a[i][j];//矩阵每行的和}cout<<b[i]<<" ";}cout<<endl;cout<<"输入精度："<<endl;cin>>s;for(k=1;; k++){for(i=0; i<n; i++){for(j=0; j<n; j++){sum=a[i][j]*x0[j]+sum;}x[i]=x0[i]+((b[i]-sum)/a[i][i]);//雅克比迭代方法计算方式c[i]=fabs(x[i]-x0[i]);//求差值的绝对值x0[i]=x[i];sum=0;}r=max(c);if(r<s)//输出迭代次数{for(i=0; i<n; i++)cout<<"x"<<i<<" = "<<x[i]<<setprecision(10)<<endl;cout<<"迭代次数："<<k<<endl;return 0;}}return 0;}GS迭代法#include<iostream>#include<math.h>#include <iomanip>#include <stdio.h>using namespace std;int n;double max(double array[100])//求最大值函数{double a=array[0];int i;for(i=1;i<n;i++){if(a<array[i])a=array[i];}return a;}main() {double s=0;double max(double array[100]);double c[100]={0.0};double x[100]={0.0};//第k+1次迭代结果，储存解数组double x0[100]={0.0};//初始向量int i,k,j;double r,sum=0;double a[100][100];double b[100]= {0.0};cout<<"请输入维数："<<endl;cin>>n;cout<<"输出a数组："<<endl;for(i=0; i<n; i++){for(j=0; j<n; j++){a[i][j]=1.0/(i+j+1);printf("%10.6lf ",a[i][j]);//矩阵中的数精确到六位}cout<<endl;}cout<<"输出b数组："<<endl;for(i=0; i<n; i++){for(j=0; j<n; j++){b[i]+=a[i][j];//矩阵每行的和}cout<<b[i]<<" ";}cout<<endl;cout<<"输入精度："<<endl;cin>>s;for(k=1;;k++){for(i=0;i<n;i++){for(j=0;j<n;j++){sum=a[i][j]*x0[j]+sum;}x[i]=x0[i]+((b[i]-sum)/a[i][i]);//gs迭代方法计算方式c[i]=fabs(x[i]-x0[i]);//求差值的绝对值x0[i]=x[i];sum=0;}r=max(c);if(r<s)//输出迭代次数{for(i=0;i<n;i++)cout<<"x"<<i<<" = "<<x[i]<<setprecision(10)<<endl;cout<<"迭代次数："<<k<<endl;return 0;}}}SOR迭代法#include<iostream>#include<math.h>#include <iomanip>#include <stdio.h>using namespace std;int n;double max(double array[100]){double a=array[0];int i;for(i=1;i<n;i++){if(a<array[i])a=array[i];}return a;}int main() {double s=0,w=0;int i,k,j;//k为迭代次数double max(double array[100]);double c[100]= {0.0}; //double x[100]= {0.0}; //第k+1次迭代结果，储存解数组double x0[100]= {0.0}; //初始向量int i,k,j;double r,sum=0;double a[100][100];double b[100]= {0.0};cout<<"请输入维数："<<endl;cin>>n;cout<<"输出a数组："<<endl;for(i=0; i<n; i++){for(j=0; j<n; j++){a[i][j]=1.0/(i+j+1);printf("%10.6lf ",a[i][j]);//cout<<a[i][j]<<" ";}cout<<endl;}cout<<"输出b数组："<<endl;for(i=0; i<n; i++){for(j=0; j<n; j++){b[i]+=a[i][j];}cout<<b[i]<<" ";}cout<<endl;cout<<"输入精度："<<endl;cin>>s;cout<<"输入松弛因子："<<endl;cin>>w;for(k=1;;k++){for(i=0;i<n;i++){for(j=0;j<n;j++){sum=a[i][j]*x0[j]+sum;}x[i]=x0[i]+(w*(b[i]-sum)/a[i][i]);//sor迭代方法的计算公式c[i]=fabs(x[i]-x0[i]);//求差值的绝对值x0[i]=x[i];sum=0;}r=max(c);if(r<s)//输出迭代次数{for(i=0;i<n;i++)cout<<"x"<<i<<" = "<<x[i]<<setprecision(10)<<endl;cout<<"迭代次数："<<k<<endl;return 0;}}}四、运行输出结果：JacobiGSSOR五、结果比较分析：说明：Hx=b，H矩阵可以由函数hi,j=1/(i+j-1)直接由程序生成，为了设定参考解，我们先设x为分量全1的向量，求出b,然后将H和b作为已知量，求x，与设定的参考解对比。