《2023年合工大计算方法实验报告A版.docx》由会员分享，可在线阅读，更多相关《2023年合工大计算方法实验报告A版.docx（45页珍藏版）》请在金锄头文库上搜索。
1、合工大计算方法实验报告A版 计算方法 试验报告 班级： 学号： 姓名： 实验一、牛顿下山法 1 实验目的 (1) 熟悉非线性方程求根简单迭代法，牛顿迭代及牛顿下山法 (2) 能编程实现简单迭代法，牛顿迭代及牛顿下山法 (3) 认识选择迭代格式的重要性 (4) 对迭代速度建立感性的认识；分析实验结果体会初值对迭代的影响 2 实验内容 （1）用牛顿下山法解方程 x 3 x 1 0 （初值为 0.6） 输入：初值，误差限，迭代最大次数，下山最大次数 输出：近似根各步下山因子 （2）设方程 f(x)=x - 3x ndash; 1=0 有三个实根 x 3 =1.8793 , x * =-0.34727
2、 ,x * * 1 2 3 或 x * 用下面六种不同计算格式，求 f(x)=0 的根 x * 1 2 3x 1 x 3 1 x = ； x = ；x = 3 3x 1 ； 2 3 x =-1.5322022 现采 x = 1 1 1 x 3 3x 1 ；x = 3 ；x = x - x 3 x 3 x 2 1 2 输入：初值，误差限，迭代最大次数 输出：近似根、实际迭代次数 3 算法基本原理 求非线性方程组的解是科学计算常遇到的问题，有很多实际背景各种算法层出不穷， 其中迭代是主流算法。只有建立有效的迭代格式，迭代数列才可以收敛于所求的根。因此设 计算法之前，对于一般迭代进行收敛性的判断是至
3、关重要的。牛顿法也叫切线法，是迭代算 法中典型方法，只要初值选取适当，在单根附近，牛顿法收敛速度很快，初值对于牛顿迭代 至关重要。当初值选取不当可以采用牛顿下山算法进行纠正。 一般迭代： x k 1 (x k ) x (x) f (x) 0 牛顿公式： x k 1 x k f (x ) k f (x ) k 牛顿下山公式： x k 1 x k f (x ) k f (x ) k 1 1 1 下山因子 1， ， ， 2 2 2 2 3 ， 下山条件 | f (x k 1 ) | f (x k )| 4 算法描述 一般迭代算法见流程图 牛顿下山算法见流程图： 图 3.1 一般迭代算法流程图 图 3
4、.2 牛顿下山算法流程图 5、代码： #include iostreamgt; #include fstreamgt; #include cmathgt; using namespace std; class srrt private: int n; double *a, *xr, *xi; public: srrt (int nn) n = nn; a = new doublen+1; /动态分配内存 xr = new doublen; xi = new doublen; void input (); /由文件读入代数方程系数 void srrt_root (); /执行牛顿下山法 void
5、 output (); /输出根到文件并显示 srrt () delete a, xr, xi; ; void srrt:input () /由文件读入代数方程系数 int i; char str120; cout n 输入文件名: ; cin gt;gt;str1; ifstream fin (str1); if (!fin) cout n 不能打开这个文件 str1 endl; exit(1); for (i=n; igt;=0; i-) fin gt;gt;ai; /读入代数方程系数 fin.close (); void srrt:srrt_root () /执行牛顿下山法 int m,
6、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 (mgt;0)amp;amp;(fabs(am)+1.0=1.0) m=m-1; if (m=0) cout n 程序工作失败！ endl; return; for (i=0; i=m; i+) ai=ai/am; for (i=0; i=m/2; i+) w=ai; ai=am-i; am-i=w; k=m; is=0; w=1.0; jt=1; while (jt=1) pq=fabs(ak); while (
7、pq1.0e-12) xrk-1=0.0; xik-1=0.0; k=k-1; if (k=1) xr0=-a1*w/a0; xi0=0.0; return; pq=fabs(ak); q=log(pq); q=q/(1.0*k); q=exp(q); p=q; w=w*p; for (i=1; i=k; i+) ai=ai/q; q=q*p; x=0.0001; x1=x; y=0.2; y1=y; dx=1.0; g=1.0e+37; l40: u=a0; v=0.0; for (i=1; i=k; i+) p=u*x1; q=v*y1; pq=(u+v)*(x1+y1); u=p-q+a
8、i; v=pq-p-q; g1=u*u+v*v; if (g1gt;=g) if (is!=0) it=1; if (it=0) is=1; dd=sqrt(dx*dx+dy*dy); if (ddgt;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; els
9、e it=1; while (it=1) t=t/1.67; it=0; x1=x-t*dx; y1=y-t*dy; if (kgt;=50) p=sqrt(x1*x1+y1*y1); q=exp(85.0/k); if (pgt;=q) it=1; if (tgt;=1.0e-03) goto l40; if (ggt;1.0e-18) it=0; if (it=0) is=1; dd=sqrt(dx*dx+dy*dy); if (ddgt;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*
10、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; xrk-1=x*w; xik-1=-y*w; k=k-1; for (i=1; i=k; i+) ai=ai-ai-1*p; ai+1=ai+1-ai-1*q; xrk-1=x*w; xik-1=y*w; k=k-1; if (
11、k=1) xr0=-a1*w/a0; xi0=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; xrk-1=x*w; xik-1=-y*w; k=k-1; for (i=1; i=k; i+) ai=ai-ai-1*p; ai+1=ai+1-ai-1*q; xrk-1=x*w; xik-1=y*w; k=k-1; if (k=1) xr0=-a1*w/a0; xi0=0.0; else u1=k*a0; v1=0.
12、0; for (i=2; i=k; i+) p=u1*x; q=v1*y; pq=(u1+v1)*(x+y); u1=p-q+(k-i+1)*ai-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 (ddgt;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
13、=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; xrk-1=x*w; xik-1=-y*w; k=k-1; for (i=1; i=k; i+) ai=ai-ai-1*p; ai+1=ai+1-ai-1*q; xrk-1=x*w; xik-1=y*w; k=k-1; if (k=1) xr0=-a1*w/a0; xi0=0.0; else dx=(u*u1+v*v1)/p; dy
14、=(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 (kgt;=50) p=sqrt(x1*x1+y1*y1); q=exp(85.0/k); if (pgt;=q) it=1; if (tgt;=1.0e-03) goto l40; if (ggt;1.0e-18) it=0; if (it=0) is=1; dd=sqrt(dx*dx+dy*dy); if (ddgt;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;
