计算方法实验报告(附代码)

实验一 牛顿下山法实验说明：求非线性方程组的解是科学计算常遇到的问题，有很多实际背景．各种算法层出不穷，其中迭代是主流算法。

只有建立有效的迭代格式，迭代数列才可以收敛于所求的根。

因此设计算法之前，对于一般迭代进行收敛性的判断是至关重要的。

牛顿法也叫切线法，是迭代算法中典型方法，只要初值选取适当，在单根附近，牛顿法收敛速度很快，初值对于牛顿迭代 至关重要。

当初值选取不当可以采用牛顿下山算法进行纠正。

牛顿下山公式：)()(1k k k k x f x f x x '-=+λ下山因子 ，，，，322121211=λ下山条件|)(||)(|1k k x f x f <+实验代码：#include<iostream> #include<iomanip> #include<cmath>using namespace std;double newton_downhill(double x0,double x1); //牛顿下山法函数，返回下山成功后的修正初值double Y; //定义下山因子Y double k; //k为下山因子Y允许的最小值double dfun(double x){return 3*x*x-1;} //dfun()计算f(x)的导数值double fun1(double x){return x*x*x-x-1;} //fun1()计算f(x)的函数值double fun2(double x) {return x-fun1(x)/dfun(x);} //fun2()计算迭代值int N; //N记录迭代次数double e; //e表示要求的精度int main(){double x0,x1;cout<<"请输入初值x0:";cin>>x0;cout<<"请输入要求的精度:";cin>>e;N=1;if(dfun(x0)==0){cout<<"f'(x0)=0,无法进行牛顿迭代!"<<endl;}x1=fun2(x0);cout<<"x0"<<setw(18)<<"x1"<<setw(18)<<"e"<<setw(25)<<"f(x1)-f(x0)"<<endl;cout<<setiosflags(ios::fixed)<<setprecision(6)<<x0<<" "<<x1<<" "<<fabs(x1-x0)<<" "<<fabs(fun1(x1))-fabs(fun1(x0))<<endl;if(fabs(fun1(x1))>=fabs(fun1(x0))){ //初值不满足要求时,转入牛顿下山法x1=newton_downhill(x0,x1);} //牛顿下山法结束后,转入牛顿迭代法进行计算while(fabs(x1-x0)>=e){ //当精度不满足要求时N=N+1;x0=x1;if(dfun(x0)==0){cout<<"迭代途中f'(x0)=0,无法进行牛顿迭代!"<<endl;} x1=fun2(x0);cout<<setiosflags(ios::fixed)<<setprecision(6)<<x0<<" "<<x1<<" "<<fabs(x1-x0)<<endl;}cout<<"迭代值为:"<<setiosflags(ios::fixed)<<setprecision(6)<<x1<<'\n';cout<<"迭代次数为:"<<N<<endl;return 0;}double newton_downhill(double x0,double x1){Y=1;cout<<"转入牛顿下山法,请输入下山因子允许的最小值:";cin>>k;while(fabs(fun1(x1))>=fabs(fun1(x0))){if(Y>k){Y=Y/2;}else {cout<<"下山失败!";exit(0);}x1=x0-Y*fun1(x0)/dfun(x0);}//下山成功则cout<<"下山成功!Y="<<Y<<",转入牛顿迭代法计算!"<<endl;return x1;}实验结果：图4.1G-S 迭代算法流程图实验二 高斯-塞德尔迭代法实验说明：线性方程组大致分迭代法和直接法。

只有收敛条件满足时，才可以进行迭代。

雅可比及高斯-塞德尔是最基本的两类迭代方法，最大区别是迭代过程中是否引用新值进行剩下的计算。

高斯-塞德尔迭代：⎪⎪⎭⎫ ⎝⎛--=∑∑+=-=++n i j k j ij i j k j ij i iik ix a x a b a x 1)(11)1()1(1n i ,...,2,1= 实验代码：#include <iostream> #include <fstream> #include <cmath> using namespace std; class gsdl {private:int n;double **a, *b, *x, eps;public:gsdl (int nn){int i;n = nn;a = new double*[n]; //动态分配内存空间for (i=0; i<n; i++) a[i] = new double[n];b = new double[n];x = new double[n];}void input (); //从文件读入系数矩阵A、常数向量B以及evoid a_gsdl ();void output (); //输出结果到文件并显示~gsdl (){int i;for (i=0; i<n; i++) { delete [] a[i]; }delete [] a;delete [] b, x;}};void gsdl::input () //从文件读入系数矩阵A、常数向量B以及e {int i, j;char str1[20];cout <<"\n输入文件名: ";cin >>str1;ifstream fin (str1);if (!fin){ cout <<"\n不能打开这个文件" <<str1 <<endl; exit(1); }for (i=0; i<n; i++) //读入矩阵Afor (j=0; j<n; j++) fin >>a[i][j];for (i=0; i<n; i++) fin >>b[i]; //读入常数向量Bfin >>eps; //读入efin.close ();}void gsdl::a_gsdl (){int i,j;double p,t,s,q;for (i=0; i<=n-1; i++){p=0.0; x[i]=0.0;for (j=0; j<=n-1; j++)if (i!=j) p=p+fabs(a[i][j]);if (p>=fabs(a[i][i])){cout <<"\n程序工作失败！" <<endl;return;}}p=eps+1.0;while (p>=eps){p=0.0;for (i=0; i<=n-1; i++){t=x[i]; s=0.0;for (j=0; j<=n-1; j++)if (j!=i) s=s+a[i][j]*x[j];x[i]=(b[i]-s)/a[i][i];q=fabs(x[i]-t)/(1.0+fabs(x[i]));if (q>p) p=q;}}}void gsdl::output () //输出结果到文件并显示{int i;char str2[20];cout <<"\n输出文件名: ";cin >>str2;ofstream fout (str2);if (!fout){ cout <<"\n不能打开这个文件" <<str2 <<endl; exit(1); }fout <<endl; cout <<endl;for (i=0; i<n; i++){fout <<x[i] <<" ";cout <<x[i] <<" ";}fout <<endl; cout <<endl;fout.close ();}void main () //主函数{gsdl c(3);c.input (); //从文件读入系数矩阵A、常数向量B以及ec.a_gsdl ();c.output (); //输出结果到文件并显示 }实验结果：实验三 选主元高斯消去法实验说明：消元是最简单的直接法，并且也十分有效的，列主元高斯消去法对求解一般的线性方程组都适用，同时可以用来求矩阵对应的行列式。

列主元高斯消去法：列主元;0||max ||≠=≤≤ik ni k lk a a消元)...,,1,()()()1()()()1(n k j i bm b b a m a a k kik k i k i k kjik k ij k ij +=⎩⎨⎧-=-=++ 回代)1...,,()(1)()(n i ax abx i iini j ji iji ii =-=∑+=实验代码：#include<iostream>#include <fstream>#include <cmath>using namespace std;class gaus{ private:int n;double **a, *b;public:gaus (int nn){int i;n = nn;a = new double*[n]; //动态分配内存空间for (i=0; i<n; i++) a[i] = new double[n];b = new double[n];}void input (); //从文件读入系数矩阵A以及常数向量Bvoid gauss ();void output (); //输出结果到文件并显示~gaus (){int i;for (i=0; i<n; i++) { delete [] a[i]; }delete [] a;delete [] b; }};void gaus::input () //从文件读入系数矩阵A以及常数向量B {int i, j;char str1[20];cout <<"\n输入文件名: ";cin >>str1;ifstream fin (str1);if (!fin){ cout <<"\n不能打开这个文件" <<str1 <<endl; exit(1); }for (i=0; i<n; i++) //读入矩阵Afor (j=0; j<n; j++) fin >>a[i][j];for (i=0; i<n; i++) fin >>b[i]; //读入常数向量Bfin.close ();}void gaus::gauss (){ int *js,l,k,i,j,is;double d,t;js = new int[n];l=1;for (k=0; k<=n-2; k++){ d=0.0;for (i=k;i<=n-1;i++)for (j=k;j<=n-1;j++){t=fabs(a[i][j]);if (t>d) { d=t; js[k]=j; is=i;} }if (d+1.0==1.0) l=0;else{ if (js[k]!=k)for (i=0;i<=n-1;i++){ t=a[i][k];a[i][k]=a[i][js[k]];a[i][js[k]]=t;}if (is!=k){ for (j=k;j<=n-1;j++){t=a[k][j];a[k][j]=a[is][j];a[is][j]=t;}t=b[k]; b[k]=b[is]; b[is]=t;}}if (l==0){ delete [] js;cout <<"\n系数矩阵奇异！无解." <<endl;return;}d=a[k][k];for (j=k+1;j<=n-1;j++)a[k][j]=a[k][j]/d;b[k]=b[k]/d;for (i=k+1;i<=n-1;i++){for (j=k+1;j<=n-1;j++)a[i][j]=a[i][j]-a[i][k]*a[k][j];b[i]=b[i]-a[i][k]*b[k];}}d=a[n-1][n-1];if (fabs(d)+1.0==1.0){delete [] js;cout <<"\n系数矩阵奇异！无解." <<endl;return;}b[n-1]=b[n-1]/d;for (i=n-2;i>=0;i--){t=0.0;for (j=i+1;j<=n-1;j++)t=t+a[i][j]*b[j];b[i]=b[i]-t;}js[n-1]=n-1;for (k=n-1;k>=0;k--)if (js[k]!=k){t=b[k]; b[k]=b[js[k]]; b[js[k]]=t;}delete [] js;}void gaus::output () //输出结果到文件并显示{int i;char str2[20];cout <<"\n输出文件名: ";cin >>str2;ofstream fout (str2);if (!fout){ cout <<"\n不能打开这个文件" <<str2 <<endl; exit(1); }fout <<endl; cout <<endl;for (i=0; i<n; i++){fout <<b[i] <<" ";cout <<b[i] <<" ";}fout <<endl; cout <<endl;fout.close ();}void main (){gaus c(3);c.input (); //从文件读系数矩阵A以及常数向量Bc.gauss ();c.output (); //输出结果到文件并显示}实验结果：实验四样条插值实验说明：实现样条插值。