合工大计算方法实验报告

《计算方法》试验报告班级：学号：姓名：实验一、牛顿下山法1 实验目的(1)熟悉非线性方程求根简单迭代法，牛顿迭代及牛顿下山法 (2)能编程实现简单迭代法，牛顿迭代及牛顿下山法 (3)认识选择迭代格式的重要性 (4) 对迭代速度建立感性的认识；分析实验结果体会初值对迭代的影响2 实验内容（1）用牛顿下山法解方程013=--x x （初值为0.6）输入：初值，误差限，迭代最大次数，下山最大次数输出：近似根各步下山因子（2）设方程f(x)=x 3- 3x –1=0 有三个实根 x *1=1.8793 , x *2=-0.34727 ,x *3=-1.53209现采用下面六种不同计算格式，求 f(x)=0的根 x *1 或x *2 x = 213xx +； x = 313-x ；x = 313+x ； x = 312-x ；x = x13+；x = x - ()1133123---x x x输入：初值，误差限，迭代最大次数输出：近似根、实际迭代次数 3 算法基本原理求非线性方程组的解是科学计算常遇到的问题，有很多实际背景．各种算法层出不穷，其中迭代是主流算法。

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

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

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

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

一般迭代：)(1k k x x φ=+ 0)()(=⇔=x f x x φ 牛顿公式：)()(1k k k k x f x f x x '-=+ 牛顿下山公式：)()(1k k k k x f x f x x '-=+λ图3.1一般迭代算法流程图 下山因子 ，，，，322121211=λ 下山条件|)(||)(|1k k x f x f <+4 算法描述一般迭代算法见流程图牛顿下山算法见流程图：5、代码：#include <iostream>#include <fstream>#include <cmath>using namespace std;class srrt{private:int n;double *a, *xr, *xi;图3.2牛顿下山算法流程图⇐⇐⇐⇐srrt (int nn){n = nn;a = new double[n+1]; //动态分配内存xr = new double[n];xi = new double[n];}void input (); //由文件读入代数方程系数void srrt_root (); //执行牛顿下山法void output (); //输出根到文件并显示~srrt (){ delete [] a, xr, xi; }};void srrt::input () //由文件读入代数方程系数{int i;char str1[20];cout <<"\n输入文件名: ";cin >>str1;ifstream fin (str1);if (!fin){ cout <<"\n不能打开这个文件" <<str1 <<endl; exit(1); } for (i=n; i>=0; i--) fin >>a[i]; //读入代数方程系数fin.close ();}void srrt::srrt_root () //执行牛顿下山法{int m,i,jt,k,is,it;double t,x,y,x1,y1,dx,dy,p,q,w,dd,dc,c;double g,u,v,pq,g1,u1,v1;m=n;while ((m>0)&&(fabs(a[m])+1.0==1.0)) m=m-1;if (m<=0){cout <<"\n程序工作失败！" <<endl;return;}for (i=0; i<=m; i++) a[i]=a[i]/a[m];for (i=0; i<=m/2; i++){w=a[i]; a[i]=a[m-i]; a[m-i]=w;k=m; is=0; w=1.0;jt=1;while (jt==1){pq=fabs(a[k]);while (pq<1.0e-12){xr[k-1]=0.0; xi[k-1]=0.0; k=k-1;if (k==1){xr[0]=-a[1]*w/a[0]; xi[0]=0.0;return;}pq=fabs(a[k]);}q=log(pq); q=q/(1.0*k); q=exp(q);p=q; w=w*p;for (i=1; i<=k; i++){ a[i]=a[i]/q; q=q*p; }x=0.0001; x1=x; y=0.2; y1=y; dx=1.0;g=1.0e+37;l40:u=a[0]; v=0.0;for (i=1; i<=k; i++){p=u*x1; q=v*y1;pq=(u+v)*(x1+y1);u=p-q+a[i]; v=pq-p-q;}g1=u*u+v*v;if (g1>=g){if (is!=0){it=1;if (it==0){is=1;dd=sqrt(dx*dx+dy*dy);if (dd>1.0) dd=1.0;dc=6.28/(4.5*k); c=0.0;}while(1==1){c=c+dc;dx=dd*cos(c); dy=dd*sin(c);x1=x+dx; y1=y+dy;if (c<=6.29) { it=0; break; }dd=dd/1.67;if (dd<=1.0e-07) { it=1; break; }c=0.0;}if (it==0) goto l40;}else{it=1;while (it==1){t=t/1.67; it=0;x1=x-t*dx;y1=y-t*dy;if (k>=50){p=sqrt(x1*x1+y1*y1);q=exp(85.0/k);if (p>=q) it=1;}}if (t>=1.0e-03) goto l40;if (g>1.0e-18){it=0;if (it==0){is=1;dd=sqrt(dx*dx+dy*dy);if (dd>1.0) dd=1.0;dc=6.28/(4.5*k); c=0.0;}while(1==1){c=c+dc;dx=dd*cos(c); dy=dd*sin(c);x1=x+dx; y1=y+dy;if (c<=6.29) { it=0; break; }dd=dd/1.67;if (dd<=1.0e-07) { it=1; break; }c=0.0;}if (it==0) goto l40;}}if (fabs(y)<=1.0e-06){ p=-x; y=0.0; q=0.0; }else{p=-2.0*x; q=x*x+y*y;xr[k-1]=x*w;xi[k-1]=-y*w;k=k-1;}for (i=1; i<=k; i++){a[i]=a[i]-a[i-1]*p;a[i+1]=a[i+1]-a[i-1]*q;}xr[k-1]=x*w; xi[k-1]=y*w;k=k-1;if (k==1){ xr[0]=-a[1]*w/a[0]; xi[0]=0.0; }}else{g=g1; x=x1; y=y1; is=0;if (g<=1.0e-22){if (fabs(y)<=1.0e-06){ p=-x; y=0.0; q=0.0; }else{p=-2.0*x; q=x*x+y*y;xr[k-1]=x*w;xi[k-1]=-y*w;k=k-1;}for (i=1; i<=k; i++){a[i]=a[i]-a[i-1]*p;a[i+1]=a[i+1]-a[i-1]*q;}xr[k-1]=x*w; xi[k-1]=y*w;k=k-1;if (k==1){ xr[0]=-a[1]*w/a[0]; xi[0]=0.0; }}else{u1=k*a[0]; v1=0.0;for (i=2; i<=k; i++){p=u1*x; q=v1*y; pq=(u1+v1)*(x+y);u1=p-q+(k-i+1)*a[i-1];v1=pq-p-q;}p=u1*u1+v1*v1;if (p<=1.0e-20){it=0;if (it==0){is=1;dd=sqrt(dx*dx+dy*dy);if (dd>1.0) dd=1.0;dc=6.28/(4.5*k); c=0.0;}while(1==1){c=c+dc;dx=dd*cos(c); dy=dd*sin(c);x1=x+dx; y1=y+dy;if (c<=6.29) { it=0; break; }dd=dd/1.67;if (dd<=1.0e-07) { it=1; break; }c=0.0;}if (it==0) goto l40;if (fabs(y)<=1.0e-06){ p=-x; y=0.0; q=0.0; }else{p=-2.0*x; q=x*x+y*y;xr[k-1]=x*w;xi[k-1]=-y*w;k=k-1;}for (i=1; i<=k; i++){a[i]=a[i]-a[i-1]*p;a[i+1]=a[i+1]-a[i-1]*q;}xr[k-1]=x*w; xi[k-1]=y*w;k=k-1;if (k==1){ xr[0]=-a[1]*w/a[0]; xi[0]=0.0; }}else{dx=(u*u1+v*v1)/p;dy=(u1*v-v1*u)/p;t=1.0+4.0/k;it=1;while (it==1){t=t/1.67; it=0;x1=x-t*dx;y1=y-t*dy;if (k>=50){p=sqrt(x1*x1+y1*y1);q=exp(85.0/k);if (p>=q) it=1;}}if (t>=1.0e-03) goto l40;if (g>1.0e-18){it=0;if (it==0){is=1;dd=sqrt(dx*dx+dy*dy);if (dd>1.0) dd=1.0;dc=6.28/(4.5*k); c=0.0;}while(1==1){c=c+dc;dx=dd*cos(c); dy=dd*sin(c);x1=x+dx; y1=y+dy;if (c<=6.29) { it=0; break; }dd=dd/1.67;if (dd<=1.0e-07) { it=1; break; }c=0.0;}if (it==0) goto l40;}if (fabs(y)<=1.0e-06){ p=-x; y=0.0; q=0.0; }else{p=-2.0*x; q=x*x+y*y;xr[k-1]=x*w;xi[k-1]=-y*w;k=k-1;}for (i=1; i<=k; i++){a[i]=a[i]-a[i-1]*p;a[i+1]=a[i+1]-a[i-1]*q;}xr[k-1]=x*w; xi[k-1]=y*w;k=k-1;if (k==1){ xr[0]=-a[1]*w/a[0]; xi[0]=0.0; }}}}if (k==1) jt=0;else jt=1;}}void srrt::output () //输出根到文件并显示{int k;char str2[20];cout <<"\n输出文件名: ";cin >>str2;ofstream fout (str2);if (!fout){ cout <<"\n不能打开这个文件" <<str2 <<endl; exit(1); } for (k=0; k<n; k++){fout <<xr[k] <<" " <<xi[k] <<endl;cout <<xr[k] <<" +j " <<xi[k] <<endl;}fout.close ();}void main () //主函数{srrt root(6);root.input (); //由文件读入代数方程系数root.srrt_root (); //执行牛顿下山法root.output (); //输出根到文件并显示}6、输入输出输出结果如下：7、分析体会牛顿下山法作为计算方法课程中重要的知识点，在看书学习时较易大致理解其思路，但上级编写代码时却是有难度的。