数值计算方法实验课

《数值计算方法实验课》由会员分享，可在线阅读，更多相关《数值计算方法实验课（27页珍藏版）》请在装配图网上搜索。

1、实验报告课程名称：数值计算方法学期：20062007学年第二学期成绩：指导教师：李耀堂学生姓名：学生学号：实验名称：实验编号：No.实验日期：实验学时：3学院：数学与统计学院专业：数学与应用数学年级：2005级一、实验目的二、实验内容三、实验环境四. 实验方法五、实验过程1 实验步骤2 关键代码及其解释3 调试过程六、实验总结1遇到的问题及解决过程2产生的错误及原因分析3体会和收获。七：程序源代码：八教师评语：实验报告课程名称：数值计算方法学期：20062007学年第二学期成绩：指导教师：李耀堂学生姓名：学生学号：实验名称：非线性方程求根实验编号：No. 1实验日期：实验学时：3学院：数学与统

2、计学院专业：数学与应用数学年级：2005级一.试验目的：练习用数值方法求解给定的非线性方程。 二.实验内容：求解人口方程:43.5156.4 二 100e入 +(e入-1)九要求误差小于10-4。三实验环境：PC计算机，FORTRAN、C、C+、VB任选一种。 四.实验方法：牛顿法牛顿法简述：牛顿法是一种特殊的迭代法，其迭代公式为：xk+1 - xk 甘k - o12， 当数列x 收敛时，其极限值x即为方程的解。定理：给定方程f (x)二0, x g a, bD 设 f (a)f (b) 0 ；则牛顿法产生的序列x 收敛于f (x) = 0在a：b内的唯一解x。k五实验过程：1编程: 根据所用

3、算法及选用语言编出源程序2. 开机, 打开所用语言系统输入所编源程序.3. 调试程序, 修改错误直至能正确运行.4. 运行程序并输出计算结果.六实验总结: 对实验中存在的问题及解决的办法、实验结果等进行分析讨论.七、教师评语附：实验报告之例实验一 非线性方程求根姓名专业：学号：一.试验目的：练习用数值方法求解给定的非线性方程。二.实验内容：求解人口方程:156.4 二 100e入 + 435 (e入-1)九要求误差小于10-4 。三实验仪器、语言：PC计算机，C程序设计语言四. 实验方法：牛顿法牛顿法简述：牛顿法是一种特殊的迭代法，其迭代公式为xk+1 - xk 甘,k -0,1,2,， 当数

4、列(V 收敛时，其极限值x即为方程的解。定理：给定方程f (x) = 0, x e a,b1) 设 f (a) f (b) 0 ；则牛顿法产生的序列 收敛于f (x) = 0在a：b内的唯一解x。五. 实验过程:K1. 编程：用C语言编出牛顿法的源程序。2. 开机，打开C语言编译程序，键入所编程序源代码.3. 调试程序, 修改错误至能正确运行.4. 运行程序并输出计算结果.次数初值epsNXXf(xx)10.450.000001100.1528867.12060920.1530.000001100.1022540.16822630.10230.000001100.1009990.000107六

5、. 讨论：(1) 牛顿法收敛速度快，但初值不容易确定，往往由于初值取得不 当而使迭代不收敛或收敛慢，但若能保证|f(x )|f(x )(称为下山 条件)，则有可能收敛。把新的近似值看作初值的话会比原来的取得 好，有可能落入局部收敛的邻域。(2) 牛顿法要求广(兀)在x附近不为零。亦即x只能是单根,不能求重根。可用重根加速收敛法求重根。(3) 牛顿法的每一步迭代中，都要计算一次导数值，若计算广(x)比 计算函数的近似值要麻烦的多。为了避免求导数，可用差商近似代替 微商f (x ) - f (x)KK-1-x -x此时牛顿迭代法改为KK-1f (X )()(4)其真实丿、/、K f (X ) -

6、f (X ) K K-1由于人口方程来源于实际问题,代表人口增长率,值不会太大, 初值不应取得过大.否则会得到该方程的另外一个解附：程序源代码：#include#define ep le-4float f (float x) float y;y=l00*exp(x)+43.5*(exp(x)-l)/x-l56.4;return(y);float df (float x) float y; y=l00*exp(x)+43.5*( x*exp(x)-exp(x)+l)/(x*x); return(y);float root(float x) float y;if (fabs)fep)y=x;els

7、e y=root(x-f(x)/df(x);while(nN)return(y); main() float y,x0;printf(“enter x0:”); scanf(“%f”,&x0); y=root(x);printf(“%f%f,y,f(y);实验二 线性方程组的直接解法一.试验目的：练习线性方程组的直接解法。二.实验内容：求解如下三对角线方程组:_ 4-1_x1x2joo_-14-1200-14-1x3200-14- 4x4200-14x5100三实验仪器、语言：PC计算机，FORTRAN、C、C+、VB任选一种。四.实验方法：追赶法五. 实验步骤:1编程: 根据所用算法及选用语

8、言编出源程序2. 开机, 打开所用语言系统输入所编源程序.3. 调试程序, 修改错误直至能正确运行.4. 运行程序并输出计算结果.六讨论: 对实验中存在的问题及解决的办法、实验结果等进行分析 讨论.七写出实验报告：实验报告要求及书写格式见附例.附：实验报告之例实验二 线性方程组的直接解法姓名专业：学号：一.试验目的：练习线性方程组的直接解法。_ 4-1_x1x2100_-14-1200-14-1x3=200-14- 4x4200-14x5100二.实验内容：求解如下三对角线方程组:三实验仪器、语言：PC计算机，FORTRAN、C、C+、VB任选一种。 四.实验方法：追赶法追赶法简述：追赶法主要

9、用于解三对角线方程组，尤其是具有严格对角占 优的三对角线方程。用克路特分解法将三对角线矩阵bc11 abc222分解为 A=LU:a1Y2a2Y3an-1bn -1aP111cn-1bn P1 Pn -11其中Y = a , i = 2,3,，n; iia = b , P = ii ii a1a = b -y Piii i-1 P = c , i = 2,3,n -1i aia =b -y P n n n n -1求解 Ax = f 可通过求解两个三角方程组 Ly = f ,Ux = y 来实现，其中y =- ,y = (f y y )/a (i = 1,2，,n).1 a i i i i -

10、1i1x = y , x = y P x )(i = n 1, n 2, ,2,1).五实验过程:n i i i i+11. 编程：用C语言编出追赶法的源程序。2. 开机，打开C语言编译程序，键入所编程序源代码.3. 调试程序, 修改错误至能正确运行.4. 运行程序并输出计算结果.计算结果：x(1) = 46.1538x(2) = 84.6154x(3) = 92.3077x(4) = 84.6154x(5) = 46.1538六. 讨论：追赶法的中间运算没有数量级的很大变化，不会有严重的误差积 累,所以此方法是比较稳定的。但在计算过程中要求a.(i = 1,2,3,n)不 能为零，因此，并不

11、是任何三对角线方程组均可用追赶法求解。附：程序源代码：#include “math.h”#include “stdio.h” int treede(n,m,b,d) int n,m;double b,d;int k,j;float s;if(m!=(3*n-2) printf(“errn”);return(-2); for (k=0;k=0;k-)dk= dk- b3*k+1*dk+1;return(2)#include “stdio.h”#include “treede.c”main()int i;static double b13=4.0, 4.0, 4.0, 4.0, 4.0,-1.0,

12、 -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0;static double d5=100.0,200.0, 200.0, 200.0,100.0; if(treede(5,13,b,d)0)for (i=0;i=4;i+) printf(“x(%d)=%5.4n”,i+1,di);实验三 矩阵求逆一.试验目的：练习用数值方法求逆矩阵。 二.实验内容：求下面矩阵的逆矩阵：- 385_A =2- 7 4.194三实验仪器、语言：PC计算机，FORTRAN、C、C+、VB任选一种。 四.实验方法： 全主元高斯约当消去法五实验步骤：1编程: 根据所用算法及选用语言编

13、出源程序2. 开机, 打开所用语言系统输入所编源程序.3. 调试程序, 修改错误直至能正确运行.4. 运行程序并输出计算结果.六讨论: 对实验中存在的问题及解决的办法、实验结果等进行分析 讨论.七写出实验报告：实验报告要求及书写格式见附例.附：实验报告之例实验三 矩阵求逆姓名专业：学号：一.试验目的：练习用数值方法求逆矩阵。二.实验内容：求下面矩阵的逆矩阵：- 385_A =2- 7 4.194三实验仪器、语言：PC计算机，C语言。四.实验方法： 列主元高斯约当消去法列高斯约当消去法简述：高斯约当消去法是高斯消去法的另一种变种和改进。本算法与高斯消元法的区别在于：（1）不用乘数m ,改用行标准

14、化，把a （k） ikkk位置上的元素先变为1。（2）高斯消元法只是把a（k）这一列中a（k）下面行的元素消 kkkk为0,而咼斯一约当消去法则是把a（k）这一列兀素除a（k）= 1以外全部消为0。（3） kkkk高斯一约当消去法进行n次消元，把第n列也消为只剩一个元素为1,其余均为 0。因此，A丰0对于消元也是必要条件。n高斯约当消去法算法为：k = 1,2,，na J a / a (j = k, k +1,，n +1)kjkj kk i丰ja J a 一 a a (j = k, k +1,，n +1)ijij ik kjPRINTx = a五实验过程：ii+11编程: 用 C 语言编出全主

15、元高斯约当消去法的源程序2. 开机，打开C语言编译程序，键入所编程序源代码.3. 调试程序, 修改错误至能正确运行.4. 运行程序并输出计算结果.计算结果：MAT A-1 IS:0.0 2 5 50.3 9 5 70.2 8 511 0.0 6 8 10.0 5 5 30.0 9 3 60.1 0 6 40.1 4 8 90.0 2 13六. 讨论：高斯约当消去法算法具有以下优点和缺点：算法优点：不用换行，换列，不用回代，精度高。算法缺点：循环语句比较难组织，已选过主元素所在行所在列的元素不能 再被选作主元素，解向量的分量也不一定按次序排列。附：程序源代码：#include “stdlib.h

16、”#include “math.h”#include “stdio.h”int rinv(n,a)int n;double a;int*is,*js,i,j,k,l,u,v;double d ,p;is=malloc(n*sizeof(int);js=malloc(n*sizeof(int);for (k=0;k=n-1;k+)d=0.0;for (i=k;k=n-1;i+)for (j=k;kd)d=p;isk=i; jsk=j;if(d+1.0=1.0)free(is); free(is);print(“err*not invn”);return(0);if(isk!=k)for(j=0;

17、j=n-1;j+)u=k*n+j;v=isk*n+j;p=au;au=av;av=p;if(jsk!=k)for(i=0;i=n-1;i+) u=i*n+k;v=i*n+jsk;p=au;au=av;av=p;l=k*n+k;al=1.0/al;for(j=0;j=n-1;j+)if(j!=k)u=k*n+j; au=au*al;for(i=0;i=n-1;i+)if(i!=k)for(j=0;j=n-1;j+)if(j!=k) u=i*n+j;au=au-ai*n+k*ak*n+j;for(i=0;i=0;k-) if(isk!=k)for(i=0;j=n-1;i+) u=i*n+k;v=

18、i*n+isk;p=au;au=av;av=p;free(is); free(is);return(1);#include “stdio.h#include “rinv.c” main() int i,j;static double a33=-3,8,5, 2,-7,4,1,9,-6;double b33;for (i=0;i=2;i+)for (j=0;j=2;j+) bij =aij;i=rinv(3,a)if(i!=0)printf(“MAT A IS:n”);for (i=0;i=2;i+) for (j=0;j=2;j+) printf(“%5.4f”,bij); printf(“n

19、”);printf(“n”); printf(“MAT A-1 IS:n”);for (i=0;i=2;i+) for (j=0;j=2;j+) printf(“%5.4f”,bij); printf(“n”);printf(“n”);实验四 线性方程组的迭代解法一.试验目的：练习线性方程组的迭代解法。二.实验内容：.分别用雅可比(Jacobi)迭代法和高斯一塞德尔 (GaussSeidel)迭代法求解下列线性方程组，且比较收敛速度, 要求当|xk+1 - Xk ” 0.01是迭代终止。210x X 二 91 21) v x +10 x 2 x = 71 2 34 x + x = 623_ 4

20、10100 -x1x0 _14101052014001x3x010041064010141x25001014x66三实验仪器、语言：PC计算机，FORTRAN、C、C+、VB任选一种。四实验方法：雅可比(Jacobi)迭代法、高斯一塞德尔(GaussSeidel) 迭代法.1编程: 根据所用算法及选用语言编出源程序2. 开机, 打开所用语言系统输入所编源程序.3. 调试程序, 修改错误直至能正确运行.4. 运行程序并输出计算结果.六讨论: 对实验中存在的问题及解决的办法、实验结果等进行分析 讨论.七写出实验报告：实验报告要求及书写格式见附例.附：实验报告之例实验四 线性方程组的迭代解法姓名专业

21、：学号：一.试验目的：练习线性方程组的迭代解法。二.实验内容：.分别用雅可比(Jacobi)迭代法和高斯一塞德尔 (GaussSeidel)迭代法求解下列线性方程组，且比较收敛速度, 要求当 |x k +1 X k | 0.01是迭代终止。210x X 二 91 21) v x +10 x 2 x = 7三.实验仪器、语言： PC 计算机，1 2 3_ 410100 -x10 _-1410101x5201400-1x03-100410x6401014-1x-25001014x664 x + x = 623C 语言。四实验方法：雅可比(Jacobi)迭代法、高斯一塞德尔(GaussSeidel)

22、 迭代法.1) 雅可比方法：用迭代法解线性方程组Ax二b，设A丄)非奇异，且对角线元素 ija丰0 (i二1,2,n),把A分裂成三个矩阵之和A=L+D+U，其中0a0 aa a1112131na0a0aa2122232naa0,Da,U0 a3132333naa a0a0n1n2n,n-1nnL=a x(k) ) ij jii(i = 12，n)x(k +1)二丄(b hi a i丿TJ知ii雅可比迭代的矩阵形式为：x(k)二 B x(k) + f 其中 B 二一D-1(L + U),f 二 D-1b2)高斯塞德尔迭代法：高斯塞德尔迭代法的分量形式为：x (k+1)=/ a x 3) 一乙

23、a x( k)(i = 1,2, n).iij =1j=i高斯塞德尔迭代法的矩阵形式为：x( k+1)= Gx(k)+ f， 其中 G=-( D + L)-1U, f = (D + L) - ib 0 1 0 1五.实验过程:1. 编程：用C语言编出追赶法的源程序。2. 开机，打开C语言编译程序，键入所编程序源代码.3. 编译、运行上述程序源代码，依提示输入方程组 1）的系数 矩阵及右端向量，记录运算结果。4. 对于方程组2），将上述程序源代码中有注释的代码中的3改 为 6，再编译、运行，依提示输入方程组2）的系数矩阵及右 端向量，记录运算结果。计算结果：1) The approximate

24、solution by Jacobin Method isx=0.9975700.999190x=0.9998470.996760after 5 iterations.The approximate solution by G-S Method is0.999830 0.999939after 4 iterations.2) The approximate solution by Jacobin Method is0.9980061.9943590.9980061.9960110.9971801.996011after 12 iterations.The approximate solutio

25、n by G-S Method is0.996554_1.9970580.998745x =1.9979200.9982251.999242after 7 iterations.六. 讨论：由于高斯塞德尔方法与雅可比法迭代阵不同，雅可比方法收敛并不能保 证高斯塞德尔方法收敛，反之也如此。但运算结果表明当二者均收敛时，高 斯塞德尔方法比雅可比方法收敛速度快。另一方面，高斯塞德尔方法算出一 个新的分量就把先一次的分量覆盖，只需一组工作单元，算法更加简单。附：程序源代码：#define N 3 /*依具体情况，3 可改为相应方程组的阶数*/#includemath.hvoid jacobi2(fio

26、at mN N, float fN, float eps, long imax)double a,b,x,eps; int i, j, tail=0;long k;float err, sum xN, xkN;for(i=0;iN;i+) xki=0;k=1;while (k=imax)err=0;for(i=0;iN;i+0)sum=0;for (j=0;j=N;j+)if(j!=i) sum+=mi j*xkj;xi=(fi-sum)/mi i;err+=(xi-xki)*(xi-xki);for(i=0;iN;i+) xki=xi;if (sqrt (err)=eps)printf(“T

27、he approximate solution by Jacoban method isn”);for(i=0;iN;i+0) printf(“%fn”,xi);printf(“after %1d iterations. n”, k);tail=1;break;k+;if (!tail) printf(“nMaximum number of iteration exceeded. n”);void gauss_siedel2(float mN N, float fN, float eps, long imax) int i, j, tail=0;long k;float err, sum, x

28、tmp, xN;for(i=0;iN;i+) xi=0;k=1;while (k=imax)err=0;for(i=0;iN;i+)xtmp=xi;sum=0;for (j=0;j=N;j+)if (j!=i) sum+=mi j*xj;xi=(fi-sum)/mi i; err+=(xi-xtmp)*(xi-xtmp);if (sqrt (err)=eps)printf(“The approximate solution by G-S method isn”); for(i=0;iN;i+0) printf(“%fn”,xi); printf(“after %1d iterations. n

29、”, k); tail=1; break;k+;if (!tail) printf(“nMaximum number of iteration exceeded. n”); main()int I, j;float aN N, bN, trans;printf(“Please input the coefficiency matrix. n”); for(i=0;iN;i+)for (j=0;j=N;j+) scanf(“%f”, &trans);ai j=trans;printf(“Please input the ordinary vector. n”); for(i=0;iv,就把a作为

30、主元素进行平面旋转变换，若a iv, ij 1ijij 1就让 a 过关一遍扫描结束后继续进行第二遍扫描(原因是原来满足 ija v的过关元素，经过其他变换还有可能过不了关)，直到扫到某 ij 1一遍时，全都过关，再设第二道关口，v =，重复前面的过程，又 2n再设v ,v ,等关口，直到v ()v为止，其中p是精度要求.34t n 0五.实验过程:1. 编程：用C语言编出雅可比(Jacobi)过关法的源程序。2. 开机，打开C语言编译程序，键入所编程序源代码.3. 调试程序, 修改错误至能正确运行.4. 运行程序并输出计算结果.计算结果：九=2.5365,对应的特征向量为x = (0.531

31、5,0.4615,0.7103)t ;11九=1.4801,对应的特征向量为x = (0.4443,0.5621,-0.6976)t ;22九=-0.0166,对应的特征向量为x = (0.7212,-0.6863,-0.0937)t 33六：讨论：雅可比过关法是雅可比方法的改进，是为消除雅可比方法因要 选主元素而花费很多机器时间而设计的，从上述计算可看出雅可比 过关法确有此优点。附：程序源代码：新建文件jcbj.c 如下：#include”math.h”void jcbj(n,a,v,eps)int n;double a,v,eps;int i, j,p,q, u,w,t,s;double

32、ff,fm,cn,sn,omega,x,y,d; for(i=0;in-1;i+)vi*n+i=1.0; for (j=0;j=n-1;j+) if(i!=j)vi*n+j=0.0;ff=0.0;for(i=0;in-1;i+) for (j=0;j=n-1;j+)d=ai*n+j;ff=ff+d*d; ff=sqrt(2.0*ff); loop0: ff=ff(1.0*n);loop1: for(i=1;in-1;i+) for (j=0;jff)p=i; q=j; goto loop; if(ffeps)return; goto loop0: loop:u=p*n+1;w=p*n+p;s=

33、q*n+q;x=-au;y=(as-aw)/2.0;omega=x/sqrt(x*y+y* y);if(y0.0)omega=-omega; sn=1.0+sqrt(1.0+-omega*omega); sn=omega/sqrt(2.0*sn);fm=aw; aw=fm*cn*cn+as*sn*sn+au*omega; as=fm*sn*sn+as*cn*cn-au*omega; au=0.0;at=0.0;for(j=0;j=n-1;j+) if(j!=p)&(j!=q)u=p*n+j;w=q*n+j; fm=au;au=fm*cn+aw*sn; aw=-fm*sn+sn+aw*cn; f

34、or(i=0;in-1;i+) if(i!=p)&(i!=q) u=i*n+p;w=i*n+q;fm=au; au=fm*cn+aw*sn; aw=-fm*sn+sn+aw*cn;goto loop1;建立主函数文件“jcbjO.c”如下：#include”stdio.h”#include”jcbj.c”main( )int i;j;double eps;static double a33= 1.0, 1.0,0.5,1.0, 1.0, 0.25, 0.5, 0.25, 2.0 static double x7 ,b6=0.0,5.0,0.0,6.0,-2.0,6.0; eps=0.00000

35、1;jcbj(3,a,v,eps);for(i=0;i=2;i+)printf(“% 5,4fn”,ai i); printf(“nn”);for(i=0;i=2;j+)printf(“%5.4f”,vij); printf(“n”);实验六 函数的插值一.试验目的：练习用数值方法计算进行函数插值。二.实验内容：已知正弦函数表xk0.50.70.91.11.31.51.71.9sin xk0.47940.64420.78330.89120.96360.99750.99170.9463试分别计算sin（0.74）与sin（1.6）的近似值（用线性插值和抛物线插值）， 并估计它的误差。三实验仪器、

36、语言：PC计算机，FORTRAN、C、C+、VB任选一种。四.实验方法：线性插值和抛物线插值法。五实验步骤：1编程: 根据所用算法及选用语言编出源程序2. 开机, 打开所用语言系统输入所编源程序.3. 调试程序, 修改错误直至能正确运行.4. 运行程序并输出计算结果.六讨论: 对实验中存在的问题及解决的办法、实验结果等进行分析 讨论.七写出实验报告：实验报告要求及书写格式见附例.附：实验报告之例实验六 函数的插值姓名专业：学号：一.试验目的：练习用数值方法计算进行函数插值。 二.实验内容：已知正弦函数表xk0.50.70.91.11.31.51.71.9sin xk0.47940.64420.

37、78330.89120.96360.99750.99170.9463试分别计算sin（0.74）与sin（1.6）的近似值（用线性插值和抛物线插值）， 并估计它的误差。三实验仪器、语言：PC计算机，FORTRAN、C、C+、VB任选一种。四. 实验方法：线性插值和抛物线插值法。 方法简述：五. 实验步骤:1编程: 根据所用算法及选用语言编出源程序2. 开机, 打开所用语言系统输入所编源程序.3. 调试程序, 修改错误直至能正确运行.4. 运行程序并输出计算结果. 计算结果：六讨论:附：程序源代码实验七 数值积分一.试验目的：练习用数值方法计算定积分。二.实验内容： 用复化梯形求积公式和复化 S

38、impson 求积公式计算J - dx,(n = 8)4 + x 2并估计误差。J dx ,(n = 5),1 + sinx0三实验仪器、语言：PC计算机，FORTRAN、C、C+、VB任选一种。四.实验方法：复化梯形求积公式和复化 Simpson 求积公式。五实验步骤：1编程: 根据所用算法及选用语言编出源程序2. 开机, 打开所用语言系统输入所编源程序.3. 调试程序, 修改错误直至能正确运行.4. 运行程序并输出计算结果.六讨论: 对实验中存在的问题及解决的办法、实验结果等进行分析 讨论.七写出实验报告：实验报告要求及书写格式见附例.附：实验报告之例实验七 数值积分姓名专业：学号：一.试验目的：练习用数值方法计算定积分。二.实验内容： 用复化梯形求积公式和复化 Simpson 求积公式计算J - dx,(n = 8)4 + x 2并估计误差。J dx ,(n = 5),1 + sinx0三实验仪器、语言：PC计算机，FORTRAN、C、C+、VB任选一种。四.实验方法：复化梯形求积公式和复化 Simpson 求积公式。 方法简述：五.实验步骤：1编程: 根据所用算法及选用语言编出源程序2. 开机, 打开所用语言系统输入所编源程序.3. 调试程序, 修改错误直至能正确运行.4. 运行程序并输出计算结果. 计算结果：六讨论: 附：程序源代码：