一、问题描述
用SOR法求解方程组Ax=b, 其中
要求程序中不存系数矩阵A,分别对不同的阶数取w=1.1, 1.2, ...,1.9进行迭代,记录近似解x(k)达到||x(k)-x(k-1)||<10-6时所用的迭代次数k,观察松弛因子对收敛速度的影响,并观察当w <= 0或 w>=2会有什么影响?
二、计算结果与分析
(1) 在收敛标准||x(k)-x(k-1)||<e=10-6且最大迭代次数200前提下,分别取阶数n=5、20、35、50、65、80,取w=1.1, 1.2, ...,1.9进行SOR法迭代求解,计算结果见图1:
图1不同阶数下松弛因子w对收敛速度的影响(w=1.1, 1.2, ...,1.9)
从图1的计算结果可以看出,在w=1.1, 1.2, ...,1.9范围内,不同阶数时SOR法的收敛速度均随收敛松弛因子w的增加而减慢,但在该参数条件下SOR法均能达到收敛要求。
(2) 为观察w£0或w³2时w对收敛速度的影响,在收敛标准||x(k)-x(k-1)||<e=10-6且最大迭代次数250前提下,分别取阶数n=5、20、35、50、65、80,分别取w=-0.5, -0.4, ...,0和w=2, 2.1, ...,2.5进行SOR法迭代求解,计算结果见图5、图6:
图2 不同阶数下松弛因子w对收敛速度的影响(w=-0.5, -0.4, ...,0)
图3 不同阶数下松弛因子w对收敛速度的影响(w=2, 2.1, ...,2.5)
从图2和图3可看出,当松弛因子w£0或w³2时,不同阶数下SOR法均迭代到了最大收敛次数,这实际上意味着该参数条件下SOR法是发散的(w=0时迭代次数为1是因为此时实际上迭代公式的残量为0,迭代无意义)。
图4给出了在w<0和w>2时SOR法求解的两个特例,显然,迭代到250次时解已趋于无限,这正式迭代发散的表现。
图4 w<0和w>2时SOR法求解的两个特例
综上所述,对于题目给出的线性方程组来说,当松弛因子w=1.1, 1.2, ...,1.9时,SOR法迭代求解是收敛的,而当w£0或w³2时SOR法是发散的。这实际上与教材中的定理很好的符合,即0<w<2是SOR法收敛的必要条件。
三、源代码
#include "stdafx.h" #include <math.h> #include<iostream> using namespace std; #include <iomanip> int SOR(int n,double w,double e=1e-6,int maxstep=100,int displaytype=0) //SOR法求解题3子函数,n:阶数,w:松弛因子,e:收敛阈值,maxstep:最大收敛步数,返回值为收敛步骤数,displaytype:是否显示计算结果,0:不显示,非零:显示 { int i,j,m; double delta; double s,t; double *px=new double [n];//分配解向量空间 for(i=0;i<n;i++)px[i]=0;//解向量赋初值 m=0;//当前跌代次数赋初值0 if (n<=2) { cout<<"阶数过小!"<<endl; return 0; } do { delta=0;//当前跌代误差赋初值0 ///////////////////第1行计算//////////////////////////// s=px[0]; t=-4*px[0]+px[1]; px[0]=(-3-t)*w/(-4)+s; delta+=(px[0]-s)*(px[0]-s);// (px[0]-s)>=0?(px[0]-s):(-px[0]+s); ///////////////////第2至n-1行计算////////////////////// for (i=1;i<=n-2;i++) { s=px[i]; t=px[i-1]-4*px[i]+px[i+1]; px[i]=(-2-t)*w/(-4)+s; delta+= (px[i]-s)*(px[i]-s);//(px[i]-s)>=0?(px[i]-s):(-px[i]+s); } ///////////////////第n行计算//////////////////////////// s=px[n-1]; t=px[n-2]-4*px[n-1]; px[n-1]=(-3-t)*w/(-4)+s; delta+=(px[n-1]-s)*(px[n-1]-s); //(px[n-1]-s)>=0?(px[n-1]-s):(-px[n-1]+s); delta=(double) sqrt((double)delta); m+=1;//跌代次数加1 } while((delta>e)&&(m<maxstep)); /////////////////////结果输出////////////////////////////// if (displaytype!=0) { cout<<endl<<"控制参数:n="<<n<<",w="<<w<<",e="<<e<<",maxstep="<<maxstep<<endl; if ((delta>e)&&(m>=maxstep))cout<<"当前参数下跌代计算发散!"<<endl; cout<<"跌代次数:"<<m<<endl<<"解向量:"; cout<<setprecision(8);//设置解向量输出精度为8位 for(i=0;i<n-1;i++)cout<<px[i]<<","; cout<<px[n-1]<<endl<<endl; } delete [] px; //if ((delta>e)&&(m>=maxstep))return -m; return m; } int _tmain(int argc, _TCHAR* argv[]) { int n=5; //阶数 double w=1.5;//松弛因子 double e=1e-6;//收敛阈值 int maxstep=250;//最大收敛步数 ////////////////不同阶数 w=1.1--1.9范围时w对收敛步数的影响/////////////////////////// for (n=5;n<85;n+=15) { cout<<endl<<"n="<<n<<endl; if (n==5) cout<<"松弛因子:"; else cout<<"松弛因子: "; for(w=1.1;w<2.0;w+=0.1)cout<<w<<" "; cout<<endl; cout<<"收敛步数:"; for(w=1.1;w<2.0;w+=0.1) { cout<<" "<<SOR(n,w, e,maxstep)<<" "; } cout<<endl; } ///////////////////////不同阶数 w<=0和w>2时w对收敛步数的影响//////////////////////// //w<=0 cout<<setprecision(3);//设置解向量输出精度为8位 for (n=5;n<85;n+=15) { cout<<endl<<"n="<<n<<endl; cout<<"松弛因子: "; for(w=-0.5;w<=-0.1;w+=0.1)cout<<w<<" "; cout<<" 0"<<endl; cout<<"收敛步数: "; for(w=-0.5;w<=0;w+=0.1) { cout<<" "<<SOR(n,w, e,maxstep)<<" "; } cout<<endl; } //w>=2 for (n=5;n<85;n+=15) { cout<<endl<<"n="<<n<<endl; cout<<"松弛因子: "; for(w=2;w<2.6;w+=0.1)cout<<w<<" "; cout<<endl; cout<<"收敛步数:"; for(w=2;w<2.6;w+=0.1) { cout<<" "<<SOR(n,w, e,maxstep)<<" "; } cout<<endl; } //特定数据点 SOR(7,-0.3, e,maxstep,1); SOR(10,2.5, e,maxstep,1); return 0; }