计算方法实验报告

计算方法实验报告 班级：信息安全 09-1 班 学号： 姓名：王义涛 合肥工业大学计算机与信息学院 2011/5/28 目录目录 实验一实验一 牛顿下山法求解非线性方程的根牛顿下山法求解非线性方程的根.3 一、实验目的.3 二、实验内容.3 三、基本原理.3 四、算法设计与实现.3 五、输入与输出.4 六、源代码.4 实验二实验二 高斯高斯赛德尔法求线性方程组赛德尔法求线性方程组.6 一、实验目的.6 二、实验内容.6 三、算法基本原理.6 四、算法设计与实现.6 五、计算用例的参考输出.8 六、源代码.8 实验三实验三 高斯消去法高斯消去法 .10 一、实验目的.10 二、实验内容.10 三、算法基本原理.10 四、算法设计与实现.10 五、计算用例的参考输出.12 六、源代码.12 实验四实验四 ROMBERG 算法算法 .14 一、实验目的 .14 二、实验内容 .14 三、算法基本原理.14 四、算法设计与实现.14 五、计算用例的参考输出.16 六、源代码.16 实验五实验五 RUNGEKUTTA 算法算法.18 一、实验目的.18 二、实验内容.18 三、算法基本原理.18 四、算法设计与实现.19 五、计算用例的参考输出.19 六、源代码.20 感悟与体会感悟与体会 .21 实验一实验一 牛顿下山法求解非线性方程的根牛顿下山法求解非线性方程的根 一、实验目的一、实验目的 (1)熟悉非线性方程求根简单迭代法，牛顿迭代及牛顿下山法 (2)能编程实现简单迭代法，牛顿迭代及牛顿下山法 (3)认识选择迭代格式的重要性 (4)对迭代速度建立感性的认识；分析实验结果体会初值对迭代的影响 二、实验内容二、实验内容 用牛顿下山法解方程 （初值为 0.6）01 3 xx 输入：初值，误差限，迭代最大次数，下山最大次数 输出：近似根各步下山因子 三、基本原理三、基本原理 牛顿下山公式： )( )( 1 k k kk xf xf xx 下山因子， 32 2 1 2 1 2 1 1 下山条件| )(| )(| 1kk xfxf 四、算法设计与实现四、算法设计与实现 流程图 关键点： while(fabs(f(x1)fabs(f(x0) cout=r n)cout迭代失败=eps) if(df(x0)=0)cout*无法迭代*m)cout迭代失败endl;break; /迭代过程 五、输入与输出五、输入与输出 x0=0.6;e=0.;m=100;n=100 输出 六、源代码六、源代码 #include #include using namespace std; 图 3.2 牛顿下山算法流程 图 double f(double x) return (x*x-1)*x-1; double df(double x) return 3*x*x-1; double newton(double x0,double eps,int n,int m) if(df(x0)=0)cout*无法迭代*fabs(f(x0) cout=r n)cout迭代失败=eps) if(df(x0)=0)cout*无法迭代*m)cout迭代失败endl;break; /迭代过程 return x1; void main() double x0,x1,e; x0=0.6;e=0.; x1=newton(x0,e,100,100); coutx=x1endl; 实验二实验二 高斯高斯赛德尔法求线性方程组赛德尔法求线性方程组 一、实验目的一、实验目的 (1)熟悉求解线性方程组的有关理论和方法； (2)能编程实现雅可比及高斯-塞德尔迭代法； (4)根据不同类型的方程组，选择合适的数值方法。 二、实验内容二、实验内容 用 Gauss - Seidel 迭代法求解方程组 2 . 45 3 . 8210 2 . 7210 321 321 321 xxx xxx xxx 输入：系数矩阵 A，最大迭代次数 N，初始向量，误差限 e 输出：解向量 三、算法基本原理三、算法基本原理 线性方程组大致分迭代法和直接法。只有收敛条件满足时，才可以进行迭代。高斯-塞德尔 是最基本的一类迭代方法，其迭代过程中引用新值进行剩下的计算。 高斯-塞德尔迭代： n ij k jij i j k jiji ii k i xaxab a x 1 )( 1 1 )1()1( 1 ni,.,2 , 1 四、算法设计与实现四、算法设计与实现 流程图关键算法： 图 4.1G-S 迭代算法流程图 关键算法： Xi解的向量，a方程组的增广矩阵 for(k=1;k=max_k;k+) cout第k次迭代; e2=0; for(i=0;in;i+) old_x=xi; sum=0; for(j=0;jm;j+) if(i!=j) sum=sum+aij*xj; xi=(bi-sum)/aii; if(e2fabs(old_x-xi) e2=fabs(old_x-xi); coutxi ; coutendl; if(e2e1) break; 五、计算用例的参考输出五、计算用例的参考输出 六、源代码六、源代码 G_S.h #include #include using namespace std; #define MAXSIZE 50 void input(double aMAXSIZEMAXSIZE,double b,int n,int m) int i,j; cout请输入原方程的增广矩阵endl; for(i=0;iaij; cinbi; coutendl; void output(double x,int n) cout原方程组的解的向量为：endl; for(int i=0;in;i+) coutxi ; Gauss_seder.cpp #includeG_S.h int main() double aMAXSIZEMAXSIZE,bMAXSIZE,xMAXSIZE; double e1,e2,sum,old_x; int n,m,i,j,k,max_k; cout输入原方程的行数 nn; cout输入原方程的列数 mm; input(a,b,n,m); cout输入迭代初始向量xi; cout输入做大迭代次数max_k; cout输入误差上限e1; for(k=1;k=max_k;k+) cout第k次迭代; e2=0; for(i=0;in;i+) old_x=xi; sum=0; for(j=0;jm;j+) if(i!=j) sum=sum+aij*xj; xi=(bi-sum)/aii; if(e2fabs(old_x-xi) e2=fabs(old_x-xi); coutxi ; coutendl; if(e2e1) break; if(k=max_k) output(x,n); else cout超过最大迭代次数，迭代失败！endl; 实验三实验三 高斯消去法高斯消去法 一、实验目的一、实验目的 (1)熟悉求解线性方程组的有关理论和方法； (2)能编程实现列主元高斯消去法； (4)根据不同类型的方程组，选择合适的数值方法。 二、实验内容二、实验内容 用选主元高斯消去求方程组 2 . 45 3 . 8210 2 . 7210 321 321 321 xxx xxx xxx 三、算法基本原理三、算法基本原理 A. ； nn m nn n m nnn n bb b bb aa aa A.)1(. . )1(. . . . | 11 111 1 111 B.消元结果直接存储在系数矩阵中； C.当消元过程发生两行对调的情况为偶数次时，行列式值为对角线乘积， 否则为对角线乘积的相反数。 四、算法设计与实现四、算法设计与实现 列主元高斯消去法：列主元;0|max| ik nik lk aa 消元 回代 )., 1,( )()()1( )()()1( nkji bmbb amaa k kik k i k i k kjik k ij k ij ) 1.,( )( 1 )()( ni a xab x i ii n ij j i ij i i i 图 4.2 列主元的约当消去 约当消去 )., 1( / / )( )()1( )( )()1( nkj abb aaa k kk k k k k k kk k kj k kj )., 1,.1, 1.1( )()()1( )()()1( nkjnkki babb aaaa k kik k i k i k kjik k ij k ij 关键点： /选主元素 for(k=0;kn-1;k+) max=akk;max_i=k; for(i=k+1;ifabs(max) max=aik; max_i=i; if(max=0)break; if(max_i!=k)/交换两行 for(j=k;jn+1;j+) t=akj; akj=amax_ij; amax_ij=t; for(i=k+1;in;i+) aik=aik/-akk; for(j=k+1;j=0;k-) sum=0; for(j=k+1;jn;j+) sum=sum+akj*xj; xk=(akn-sum)/akk; /回代 五、计算用例的参考输出五、计算用例的参考输出 六、源代码六、源代码 Gauss.h #include #include using namespace std; #define MAX 50 void input(double aMAXMAX+1,int n) cout输入原方程组的增广矩阵endl; for(int i=0;iaij; void output(double x,int n) coutGauss 消去法得到的原方程组的解为endl; for(int k=0;kn;k+) coutxk ; GaussXiaoqu.cpp #includeGauss.h int main() double aMAXMAX+1,xMAX,sum,max,t; int n,i,j,k,max_i; cout输入原方程组的阶n; input(a,n); for(k=0;kn-1;k+)/选主元素 max=akk;max_i=k; for(i=k+1;ifabs(max) max=aik; max_i=i; if(max=0)break; if(max_i!=k)/交换两行 for(j=k;jn+1;j+) t=akj; akj=amax_ij; amax_ij=t; for(i=k+1;in;i+) aik=aik/-akk; for(j=k+1;jn+1;j+) aij=aij+aik*akj; /消元 if(max=0)cout原方程组无解=0;k-) sum=0; for(j=k+1;jn;j+) sum=sum+akj*xj; xk=(akn-sum)/akk; /回代 output(x,n); coutendl; return 0; 实验四实验四 Romberg 算法算法 一、一、实验目的实验目的 (1) 熟悉梯形递推算法、龙贝格算法； (2) 能编程实现梯形递推算法、龙贝格算法； (3) 分析实验结果体会各种方法的精确度，建立计算机求解定积分问题的感性认识 二、二、实验内容实验内容 用龙贝格算法计算dx x x 1 0 sin 输入：积分区间，误差限 输出：序列 Tn，Sn,Cn,Rn 及积分结果（参考书本 P71 的表 2-5） 三、算法基本原理三、算法基本原理 由梯形递推公式求得梯形序列，相邻序列值作线性组合得 Simpson 序列, Simpson 序列作线性组合得柯特斯序列, 柯特斯序列作线性组合的龙贝格序列。若 |R2-R1|e,则输出 R2;否则依此类推。 梯形递推公式 1 0 2 12 )( 22 n k k n n xf hT T 加权平均公式： n nn S TT 14 4 2 n nn C SS 14 4 2 2 2 n nn R CC 14 4 3 2 3 四、算法设计与实现四、算法设计与实现 kTSCR 0R00 1R10R11 2R20R21R22 3R30R31R32R33 4R40R41R42R43 本实验，未采取书本所提供的框图算法，而是直接使用 T 表，先通过梯形递推公 式求解和，紧接着就求解 S，通过精度判断，是否继续，将所得到的 n T n T2 T，S，C，R都存储在 R中，然后通过加权平均公式： n nn S TT 14 4 2 推出 ri+1j=ri+1j-1+(ri+1j-1-rij-1)/(m-1); n nn C SS 14 4 2 2 2 n nn R CC 14 4 3 2 3 每次循环将 m=4*m,m 初值为 1。 关键点： 先用梯形递推法求解出 T 存在 r0中，再通过一循环求解 S、C、R h=b-a;n=1; tn=h*(f(a)+f(b)/2; cout.precision(6); do ri0=tn; s=0; x=a+h/2; while(xb) s=s+f(x); x=x+h; t2n=(tn/2+h*s/2); h=h/2; tn=t2n; i+; while(iMAXSIZE); 加权平均公式 for(j=1;jMAXSIZE-1;j+) m=m*4; for(i=0;iMAXSIZE;i+) ri+1j=ri+1j-1+(ri+1j-1-rij-1)/(m-1); 五、计算用例的参考输出五、计算用例的参考输出 六、源代码六、源代码 #include #include #include using namespace std; #define MAXSIZE 5 double f(double x); int main() double x,a,b,h,s,e,tn,t2n,rMAXSIZEMAXSIZE; int n,i=0,j,m=1; cout请输入积分区间边界 a 和 bab; cout输入精度e; h=b-a;n=1; tn=h*(f(a)+f(b)/2; cout.precision(6); do ri0=tn; s=0; x=a+h/2; while(xb) s=s+f(x); x=x+h; t2n=(tn/2+h*s/2); h=h/2; tn=t2n; i+; while(iMAXSIZE); for(j=1;jMAXSIZE-1;j+) m=m*4; for(i=0;iMAXSIZE;i+) ri+1j=ri+1j-1+(ri+1j-1-rij-1)/(m-1); coutTn、Sn、Cn、Rn 如下表所示endl; for(i=0;i4;i+) for(j=0;ji+1;j+) coutrij ; coutendl; for(i=0;i4;i+)coutr4i; coutendl; if(fabs(r43-r33)=e) cout积分=r43endl; else cout超出最大迭代次数！，异常退出！endl; double f(double x) if(x=0)return 1; else return sin(x)/x; 实验五实验五 RungeKutta 算法算法 一、实验目的一、实验目的 (1) 熟悉数值微分中 Rung-Kutta 方法； (2) 能编程实现 Rung-Kutta 方法； (4)明确步长对算法的影响并理解变步长的 Rung-Kutta 方法 二、实验内容二、实验内容 （1） 0 x1 10 2 y y x yy 取 h=0.1 时用 Rung-Kutta 方法求其数值解并与精确解进行比较。 输入：求解区间，初值，数值解个数 输出：数值解 三、算法基本原理三、算法基本原理 通过龙格-库塔法我们可以获得更高精度。经典龙格-库塔法即在区间xn,xn+1取四点， 并对这四点的斜率进行加权平均作为平均斜率，通过泰勒公式寻找使局部截断误差为 O(h5) （即 4 阶精度）的参数满足条件。 改进的欧拉公式： 预测 ),(),( 2 111 nnnnnn yxfyxf h yy 校正 ),( 1nnnn yxhfyy 四阶（经典）龙格-库塔公式 ),( ),( ),( ),( )22( 34 2223 1222 1 432161 hKyhxfK KyxfK KyxfK yxfK KKKKyy nn h n h n h n h n nn h nn 四、算法设计与实现四、算法设计与实现 流程图： 五、计算用例的参考输出五、计算用例的参考输出 图 5.3 经典龙格库塔算法 六、源代码六、源代码 #include #include #include using namespace std; double f(double x, double y) return (y-2*x/y); void R_K(double x0,double y0,double h,int N) cout.precision(7); double x1,y1,K1,K2,K3,K4; int i=1; while(iN+1) x1=x0+h; K1=f(x0,y0); K2=f(x0+0.5*h,y0+0.5*h*K1); K3=f(x0+0.5*h,y0+0.5*h*K2); K4=f(x0+h,y0+h*K3); y1=y0+h*(K1+2*K2+2*K3+K4)/6; coutxi=;coutx1 yi=;couty1endl; i=i+1; x0=x1; y0=y1; int main() double a,b,h; double y0=1; int N; cout请输入微分范围（a,b）ab; cout请输入数据组数 NN; h=(b-a)/N; R_K(a,y0,h,N); return 0; 感悟与体会感悟与体会 学完了计算方法（数值分析） ，也做完了所有的实验，感受颇深。 这些实验，按照指导书基本都能实现，但是龙贝格算法的话，按照流程图比较复杂， 在梯形递推公式求解 T 的时候就循环一次，其自身还有好几个循环，很容易乱。但是按照 T 表就没那么复杂，先以一个循环求解出 T，存储起来，再以一个循环用加权平均公式求解 出其他几个，这样便大大减小了时间复杂度。 我想在做实验的时候，我们不能抱有侥幸的心理，要扎扎实实的做好每一个细节，如 此才能有效率的完成实验，否则，每次都要回头检查代码，浪费大量的时间。当然，做好 预习工作也是很重要的，比如说事先就写好了算法或是流程图，在写代码时就不会乱了。 做实验，关键还是老师的指导，在汪老师的指导下，注意到了以前写代码时没有注意 到的一些细节问题，有的时候细节是注定成败的。我觉得在以后也应该像这样，不清楚的 直接问老师，从老师的回答中要比我们自己看要强的多。总的来说，在完成这次试验后， 我的编程基本功又增强了一些。就像汪老师在上课所强调的那样，平时多练练这些题，编 程能力肯定是会增强的。 还有就是，做实验的时候切记不能急躁，往往急躁的时候就很容