数值计算方法

1、1.题目造倒数表，并例求18的倒数。(精度为0.0005)2算法原理2.1牛顿迭代法牛顿迭代法是通过非线性方程线性化得到迭代序列的一种方法。对于非线性方程f X()= 0，若已知根X*的一个近似值Xk，将f(X)在Xk处 展成一阶泰勒公式后忽略高次项可得：f (x) f X( k ) + f (Xk )(X - Xk )右端是直线方程，用这个直线方程来近似非线性方程f (X)。将非线性方程fx( ) = 0 的根 X*代入 f x( *) = 0，即卩f x( k ) + f (Xk )(X* - Xk ) 0-f (xk )解出X心 f (Xk )将右端取为Xk+1，贝U Xk+1是比Xk更

2、接近于X*的近似值，即f欲)Xk+ 1 Xk -f (Xk )这就是牛顿迭代公式，相应的迭代函数是f (X)?(x) = X -f (X)2.2牛顿迭代法的应用11计-算是求CX- = 1 0的解，解出X，即得到。取c c有牛顿迭代公式CXk - 11 _Xk+1 = xk -= c c这样就失去了迭代的意义，达不到迭代的效果。1f (x) = C，CX2 - X，f (x) = CX- 1, 故重新构造方程：CX2 - X = 0 ,也-是该式的解。故取f (x)Cf (x) = 2CX - 1，则有牛顿迭代公式Xk+1 = Xk - CXk2 - Xk = CXk2, k = 0,1,.2

3、CXk - 12Ck - 1丄 11的值一之间，取初值xo = 0.1。20 1013流程图开始读入Xo,，N錯東384输出结果fr OUt pLFtl txt -记孝丰划it旧離査看(V) SStofH)*捏1数表*水奪* k-0k（0）0.10000k=lx（l）=0.06923k二2x（2）=0.05781k-3x（3）=0.05564k二4x（4）=0.055641/18.000000=0, 056Jk5. 结果分析当 k= 3 时，得 5 位有效数字 0.05 564。此时，X3 - X4 = 0.00 000 X*时必收敛，但是当X0 0)时迭代结果发散，?较小尚不确定。6. 心得

4、体会起初对题目的理解并不是很透彻，另外对构建牛顿迭代公式理论依据不是特 别充分，比如说为什么在原有直接得到的式子两边各乘一个X，只是试出来的。在编程方面不够成熟。当然也加深了对牛顿迭代法的理解和应用的具体实现。实验二例3-41.题目用列主元消去法求解方程组?12xi - 3x2 + 3x3 =15? ?- 18xi - 3x2 - X3 = - 15?1 + X2 + X3 = 6?并求出系数矩阵 A的行列式的值det Ao2.算法原理2.1 顺序高斯消去法 顺序高斯消去法是利用线性方程组初等变换中的一种变换，即用一个不为零 的数乘一个方程后加至另一个方程，使方程组变成同解的上三角方程组，然后

5、再 自下而上对上三角方程组求解。这样，顺序高斯消去法可分成“消去”和“回 代” 两个过程。在用顺序高斯消去法时，在消元之前检查方程组的系数矩阵的顺序主子式， 当阶数较高时是很难做到的。若线性方程组的系数具有某种性质时，如常遇到的 对角占优方程组，自然能够用高斯消去法求解。2.2 列选主元消去法线性方程组只要系数矩阵非奇异，就存在惟一解，但是按顺序消元过程中可 能出现主元素 akk( )k = 0，这时尽管系数矩阵非奇异，消元过程无法再进行，或者即使akk( )k丸，但如果其绝对值很小，用它作除数也会导致其他元素的数量级急剧增大和使舍入误差扩大，将严重影响计算的精度。为避免在校园过程确定乘数时的

6、所用除数是零或绝对值小的数，即零主元或 小主元，在每一次消元之前，要增加一个选主元的过程，将绝对值大的元素交换 到主对角线的位置上来。列选主元是当高斯消元到第 k 步时，从 k 列的 akk 以下(包括 akk )的各元素 中选出绝对值 大的，然后通过行交换将其交换到akk 的位置上。交换系数矩阵中的两行(包括常数项) ，只相当于两个方程的位置交换了，因此，列选主元不影 响求解的结果。列选主元消去法常用来求行列式。设有矩阵?a11 L a1n ?A =?M M ?an1 L ann ?用列主元消去法将其化为上三角形矩阵，对角线上元素为01,322(2)也a33(3)，于是行列式det A =

7、- (1)ma a11(1) 22(2) Lann( )n其中 m 为所进行的行交换次数。这是实际中求行列式值的可靠方法。3流程图/读入数据心Ai, j = l,工 bnr -nr1 1 1-v Xn h* 1/= !J2 w1 *a愉岀b|jby“ J人 /结束Jisim斯消去法算法框图4输出结果outpi_jt2.txt -记事本= I 回菊琢的梅式的SBlV)x (1) = 1. 0000 x(2)=2. 0000 x(3)=3.0000 detA=- G6. 00005.结果分析采用计算机运算在计算大数据时有明显的优点，另外也需要考虑到存储。高斯消去法的使用条件是akk()k工0,k

8、= 1,2,L, n,而列选主元法可以保证这一条件。并且可以避免在消元过程确定乘数时所用除数是绝对值小的数，相对全选 主元的运算量小，一般也可以满足精度要求。6心得体会此次上机不仅需要对原理了解透彻，而且要求的编程能力较强。在定义和思 路上没问题，只是在编程软件的使用上遇到些不稳定的问题，如头文件的使用。 在存储空间上得到了新的认识，另外发现了当代码多时流程框图的好处。编程是 一件很需要耐心的事，自己还有很大进步空间。实验三例3-101.题目用杜里特尔分解法求矩阵A的逆矩阵A-1?11-1?A = ?12 -2?-211 ?2算法原理2.1杜里特尔分解法 设线性方程组Ax b=，对系数矩阵A进

9、行除不交换两行 位置得初等行变换相当于用初等矩阵 M1左乘A,在对方程组第一次消元后，A(1) 和b分别化为A(2)和?M1 A(1) = A(2)?b(2) ，即?1?-m21 1?其中M1 = ?-m311?MO?-mn 11?M1b(1) = b(2)第 k 次消元时， A( )k 和 b( )k 分别化为 A(k+ 1) 和 b(k+ 1) ，即?M Ak ( )k = A(k+ 1) ?M b( )k = b(k+ 1)?1其中?O?1M1 = ?- mk+ 1 ,k?MO?- mnk?1?消元过程是对 k=1 n- 1进行的，因此有?Mn- 1LM M A21 (1) = A( )

10、n?M - LM M b2 (1) = b( )nn 11将上三角形矩阵A()n记为U，于是有A = M M1-12- 1LMn- 11U = LU?1?m211?mm1?其中L = M M1- 12- 1LMn- - 11 = ?3132? ?MMm43 O?MMMO1?mn1mn2mn3 L?为单位下三角形矩形。这样高斯消去法的实质是将系数矩阵 A 分解为两个三角形矩阵 L 和 U 相 乘，即 A=LU在上述矩阵描述中遇到了下三角形矩阵运算。主对角线以上元素全为零的 方阵称为下三角形矩阵。下三角形矩阵的乘积仍是下三角形矩阵。若下三角形 矩阵可逆，其逆矩阵仍是下三角形矩阵，而且下三角形矩阵的

11、乘积和逆矩阵很 容易求得。把 A 分解成一个单位下三角阵和一个上三角阵 U 的乘积成为杜里特尔分 解。这种分解是惟一的。2.2 高斯 -约当法高斯消去法有消元和回代两个过程，当对消元过程稍加改变便可以使方程组 化为对角形方程组Dx b= 的形式，其中矩阵 D 为对角形矩 阵，即?a11(1)?(2)?D = ?a22? O ?ann( )n ?当高斯-约当消去法消元的每一步都先用主元去除其所在行的各元素(包括常 数项)时，方程组便可化成?1?x1 ? ?b1( )n ? 1 O?x?M? 2 ?= ?bM2( )n ? 1?xn ?( )n ? bn ?这是等号右端即为方程组的解。高斯 -约当

12、消去法每一步都用主元去除其所在行 的各元素(包括常数项) ，这个个过程成为归一化，这时方程组的系数阵转化为 单位阵。为减小误差，高斯 -约当消去法还常用列选主元技术。3 流程图Inf1=%r + =/1+ =14输出结果5.结果分析采用杜里特尔分解法求解方程组时，由于把对系数矩阵的计算和对右端项的 计算分开，这使计算线性方程组系非常方便。只需进行一次矩形三角分解，然后 再解多个三个方程组，且多解一个方程组仅需要增加大约 n2次乘除法运算。采用高斯约当法仅需要进行消元归一，而不需要回代，为编程实现提供了便 利。6心得体会步骤很重要，审题-确定算法-解题步骤-流程图-程序-代入简单值进行验 证。在

13、编程时先在代码输入区打好框架，并且尽量在每一命令后注释，方便检查 错误和日后复习。定义和变量存储很灵活，如我把单位向量直接赋给了A矩阵变量中，还有根据 终的目的直接简化计算。另外赋值前，确定存储空间并且要 定义初值为零。实验四例4-61.题目已知f（X）的观测数据X1234f( )x0-5-63构造插值多项式2算法原理首先构造基函数I Xk ()= n=n0xxk-Xfi,可以证明基函数满足下列条件：i k工?)i斗I Xk( i ) = ? ,?l i = k对于给定n+1个节点，n次拉格朗日插值多项式由下式给出:nL x()=习 x yk ()kk=0其中I Xk ( )= n=noXXk

14、 - XXii由于I Xk()疋个关于x的n次多项式，所以L x()为关于x的不咼于n次的代数多项式。当x = Xi时，L x( i ) = yi ,满足插值条件。3流程图开始输入（x时）ii = 0,1, l_,n0?y0_ k?x x-t? txk xjj = 0, ?水1,k1, + ?ny +ty? y1F输出y k+l-k4输出结果output4.ibd:-迅事耳1 = 1 Zb-I交傾F)號棹衣宜看朗甚啟z=2. 5000000处的函数值为：尸-6. 37500005.结果分析由于所知的拉格朗日计算机算法只能实现计算某一特定值的近似函数值，而 不知如何导出表达式，故例求x=2.5处

15、的函数值以说明表达式以得出，只是在计算机程序中。并且也能达到拉格朗日插值法使用的目的。6心得体会编程不够细心，程序没问题，却因为不知道是输入文件错了而检查了好长时 间。但同时也加深了对拉格朗日基函数性质的认识和理解。实验五习题5-21.题目给出平面函数 z x y(,) = ax +by + c的数据i12345X0.10.20.40.60.9yi0.20.30.50.70.8Zi0.580.630.730.830.92按 小二乘原理确定 a b c,，2算法原理2.1小二乘原理设已知某物理过程y = f x()的一组观测数据(Xi , f x( i ), i = 1,2=m要求在某特定函数类

16、 ()x中寻找一个函数?( )x作为y = f x()的近似函数，使得二者在Xi 上的误差或称残差S i = ?(x) - f x( i ) , i =1,2,Lm按某种度量标准为 小,这就是拟合问题。要求残差S i按某种度量标准为 小,即要求由残差S i构成的残差向量S = S 0,S S丄m T的某种范数S为小。例如，要求II sll,或S 即mm11 S i=x x = ?(为)-f x( i=)IIS = max S i = max ?( )x - f x()ii为小,这本来都是很自然的，可是计算不太方便。所以通常要求:或者Il sii 2 =?xm si2 ?2 二?xm ?()x-

17、 f x()i 2 ?i= 0? i =0II sll 22= Xi=0X2 = ?(x ) - f x( i ) 2i=0为 小。这种要求误差平方和 小的拟合称为曲线拟合的 小二乘法。2.2多变量数据拟合对于给定的一组数据(xi, yi ) , i=1, 2，m，寻求做n次多项式y = xa x/k= 0使性能指标J a a 0,1 丄an ) = 乂(yi - 士a Xk ik )2 为小。i=1k=0由于性能指标J可以被看做关于ak, k=0, 1,，n的多元函数，故上述拟 合多项式的构造问题可转化为多元函数的极值问题。令?J=0从而有正则方程组?m刀址Xi2?Xn?ao ?X y ?n

18、+ 1?a13l?EXy ?无刀刀炉LXiL口? M =? M?MMMM? ?ak ?XinXin+ 1 Xin+2 L C2n?对多变量（或称多元）线性模型y* = ao + a xi 1 +a X22 + L+ a xn n进行了 m次观测劝*= ao + a xi 11 +a X221 + L+a xni n?y2* = a0 +a x1 12 +a x222 +L+a xn2 n?M?y?n* = a0 +a x1 1n +a x22n +L+a xnn n这个称为回归方程组，写成矩阵形式y=Aa?y1* ?1x11x21 L xn1 ?a0 ?x22?* ? ?其中?y=?2 ?,L

19、 xn2 ?, a =?a1A=?1 x12y?。?M ? ?MMM M M ?M ?* ?1?ym ?x1mx2m L xnm ?an ?当m n时，要确定一组a ii, = 0,1,L,n,使之精确地满足 m个方程，这是超定方程组的问题，只能在 小平方误差的基础上确定 a i。定义残差向量 = S3 S 1, 2丄m T ,则S = y-A a其中 y = y y1,2,L, ym T 代替输出向量。取性能指标J =S TS使之小,以此确定出a。由J =S TS= (y - Aa) (T y - Aa) =yyT - aTAy - yTAa +aTA AT a利用向量和矩阵的运算公式，有A

20、 AT a=AT y此即为正则方程组，当 A AT 非奇异时，可求得a =(A AT)- 1AT y3流程图5?n输入x yi,Zji = 1,2,，nXi ?yi ?-2i = 1,ATi；ATi；AT3i，nna ki ajij = 1,23naki zi ?Y kz=1X kjk? = 3丰输出参数4输出结果5结果分析曲线拟合的 小二乘法是反映所给数据点的总的趋势，并不是严格的通过每 个数据点，这样就避免了大量数据插值时需要高次多项式，同时又去掉了数据 所含的测量误差。6心得体会整理思路越来越熟练了，所以执行各个步骤也相对简单了很多。另外对原 理也加深了认识。附录：1. 造倒数表1源程序

21、：#includestdio.h#includemath.h#define N 30void main()int i; float xN,c; FILE *fp1,*fp2;fp1=fopen(input1.txt,r);fp2=fopen(output1.txt,w);fscanf(fp1,%f,&c);fscanf(fp1,%f,&x0);/ 初 值fprintf(fp2, * 倒 数 表 *n);for(i=0;iN;i+)/ 牛顿迭代法 xi+1=xi*xi*c/(2*c*xi-1);fprintf(fp2,k=%dtx(%d)=%.5fn,i,i,xi); if(fabs(xi+1-

22、 xi) 输入文件： input1.txt180.13 输出文件： “ output1.txt ”倒数表 *k=0x(0)=0.10000k=1x(1)=0.06923k=2x(2)=0.05781k=3x(3)=0.05564k=4x(4)=0.05564计算结果：1/18.000000=0.0562. 例 3-41源程序：#includeiostream #includecmath using namespace std;#define N 10 voidmain()int i,j,k,l,n; floatbN,aNN,t,d,det=1.0;/* 数据输入 */FILE *fp1,*fp

23、2; fp1=fopen(input2.txt,r); fp2=fopen(output2.txt,w); fscanf(fp1,%d,&n);for(i=0;in;i+)for(j=0;jn;j+) fscanf(fp1,%f,&aij);for(i=0;in;i+) fscanf(fp1,%f,&bi);/* 数据输入 */*高斯消去 */*消元 */*列选主元函数 */ for(k=0;kn-1;k+)/ 从第一次消元到第 N-1 次消元 d=akk;l=k for(i=k+1;ifabs(d)d=aik;i=l;if(i=n)/ 判断是否奇异，不奇异进行行交换 if(d=0)fprin

24、tf(fp2, 奇异);/如果所有行的 “首列 ”都为 0，为奇异 elseif(l!=k)/ 如果第 k 行的 “首列 ”并不是 大det=det*(-1);for(j=k;j=n;j+)/ 交换系数矩阵中的两行 t=alj;alj=akj;akj=t;t=bl;bl=bk;bk=t;/ 交换右端常向量中的两行/*列选主元函数 */for(i=k+1;in;i+)/ 第( k+1 )次消元要得到 N-(k+1) 个乘数aik=aik/akk;for(j=k+1;j=0;i-)/ 从倒数第二项开始依次回代N-1 次 t=0;for(j=i+1;jn;j+)t=t+aij*bj; bi=(bi-

25、 t)/aii;/*高斯消去 */数据输出 */ for(i=0;in;i+)/ 输出方程组的/*解fprintf(fp2,x(%d)=%.4fn,i+1,bi); for(i=0;in;i+) det=det*aii;fprintf(fp2,detA=%.4fn,det);/ 输出系数矩阵行列式的值 fclose(fp1);数据输出 */fclose(fp2);/*2 输入文件： “ input2.txt ”312 -3 3 -18 3 -11 1 115 -15 63输出文件： “ output2.txt x(1)=1.0000x(2)=2.0000x(3)=3.0000 detA=-66

26、.00003. 例 3-101源程序：#includeiostream#includecmathusing namespace std;#define N 30 voidmain()int i,j,r,k,n;float aNN=0,s;/* 数据输入 */FILE *fp1,*fp2;fp1=fopen(input3.txt,r);fp2=fopen(output3.txt,w);fscanf(fp1,%d,&n);for(i=0;in;i+)for(j=0;jn;j+)fscanf(fp1,%f,&aij);for(i=0;in;i+)/ 输入单位阵 j=i+n;aij=1J/* 数据输入

27、 */*LU分解 */ for(i=1;in;i+)/r=0 时： 1 区间不变化， 2 区间变化。ai0=ai0/a00;for(i=1;in;i+)/ 第二行列至末行列进行变化时 for(j=i;j2*n;j+)/ 第 2(r+1)-1 区间的变化行 s=0.0;for(k=0;ki;k+)/ 对应行列元素的乘积进行求和s+=aik*akj; aij=aij-s;for(j=i+1;jn;j+)/2 (r+1 ) 区间的变化列 s=0.0 for(k=0;kLUY*LU分解 */* 高斯约当法解 Ux=Y*/*for(i=1;in;i+)for(j=0;ji;j+) aij=0;/* /*

28、化a0j=a0j/a00;a00=1;/ 第一列其余行已为零for(i=1;in;i+)/(n-1) 次归一消元提取 UY 减少计算 */提取 UY 减少计算 */消元 */ for(j=0;j2*n;j+)/ 首行归 for(j=i+1;j2*n;j+)/ 第 i+1 行的各列进行归一 化aij=aij/aii; aii=1;for(r=0;ri;r+)/ 对第一行至第 i-1 行的 i 列进行置零for(j=r+2;j2*n;j+)/r 行的各列与第 r+1 行的对应列进行减法运算 arj=arj- aij*arr+1;/ 第 r+1 行归一后，乘数即为要置零处的值 arr+1=0;/ 乘

29、数用完之后置零 即可消元 */*/* 高斯约当法解 Ux=Y*/ fprintf(fp2,A 的逆矩阵为 nn);/*数据输出 */for(i=0;in;i+) fprintf(fp2,| );for(j=0;j输入文件： input3.txt31 1 -11 2 -2 -2 1 13输出文件： output3.txtA 的逆矩阵为2.0000-1.00000.00001.5000-0.50000.50002.5000-1.50000.50004. 例 4-61源程序：#includestdio.h#includemath.h#define N 50 int n,i,j,k; floatxx=

30、0.0,yy=0.0,t,xN,yN,cN,AN;main()FILE *fp1,*fp2;fp1=fopen(input4.txt,r);fp2=fopen(output4.txt,w);fscanf(fp1,%d,&n);/n=4 for(i=0;in;i+) fscanf(fp1,%f,%f,&xi,&yi);/1,0,2,-5,3,-6,4,3fscanf(fp1,%f,&xx);/ 要计算的值 for(k=0;kn;k+)/ 拉格朗日插 值t=1.0; for(j=0;jk;j+) t=(xx-xj)*t/(xk-xj);for(j=k+1;jn;j+) t=(xx-xj)*t/(xk-xj);yy=yy+t*yk; fprintf(fp2,nx=%.7f 处的函数值为:y=%.7f,xx,yy);fclose(fp1);fclose(fp2);2 输入文件： input4.txt41,02, 43, 54,32.53输出文件： output4.txtx=2.5000000 处的函数值为 :y=-6.37500005. 习题 5-21源程序：#include stdio.h#include math.h#define N 30 void main()int i,n,k,j,l;float xN,yN,zN;/ 定义输