数值计算方法上机实验报告

1、华北电力大学科技学院课程名称：数值计算方法上机实验报告牛玺童电气11k6111904010415实验造倒数表231. 题目造倒数表，并例求18的倒数。(精度为0.0005)2. 算法原理x，将f(X)在x处展2.1牛顿迭代法牛顿迭代法是通过非线性方程线性化得到迭代序列的一种方法。 对于非线性方程f (X)0,若已知根X*的一个近似值成一阶泰勒公式后忽略高次项可得:f(X) f (Xk )f '(Xk )(xXk )右端是直线方程，用这个直线方程来近似非线性方程f(X)。将非线性方程f(X)0的根X*代入f(X* )0，即解出f (Xk )f '(Xk )(XXk )0X Xkf

2、 (Xk )f '(Xk)将右端取为Xk1，则Xk1是比Xk更接近于X*的近似值，即Xk 1Xkf (Xk )f '(Xk)这就是牛顿迭代公式，相应的迭代函数是(X)Xf (X)f '(X)2.2牛顿迭代法的应用11计算_是求CX 10的解，解出X，即得到_。取f (X)CX 1， f '(X)c，cc有牛顿迭代公式Xk 1 Xk CXk 11这样就失去了迭代的意义, 故重新构造方程： CX2c c 达不到迭代的效果。1-也是该式的解。故取f (X) CX2 X，cf '(X) 2CX 1，则有牛顿迭代公式1的值在183.流程图2cx x_kkXk 1

3、xk2cx1k丄丄之间，取初值20 102cxk2C1k0.1。xo , Nt读入4. 输出结果丈锹F精«5C(OJ查看 WWH)kk k k k* * ajc*怖数表* *址*K（0）=0. 100003tl）=o. 06923X（2）=0. 05781x3）=0. 05564x4）=0. 05564计篦结果=1/18. 0000000. 0565. 结果分析当k 3时，得5位有效数字0.05 564此时，X3 X4 0.00 000 0.0 005,故取 X* X30.05 564 0.056。此种迭代格式仍存在一定的缺陷，经实验后发现当初值 X0 X时必收敛，但是当X0 X*(

4、0)时迭代结果发散， 较小尚不确定。6. 心得体会起初对题目的理解并不是很透彻，另外对构建牛顿迭代公式理论依据不是特 别充分，比如说为什么在原有直接得到的式子两边各乘一个X,只是试出来的。在编程方面不够成熟。当然也加深了对牛顿迭代法的理解和应用的具体实现。实验例 3-41.题目用列主元消去法求解方程组12x13x23x31518x13x2x315x1x2x36并求出系数矩阵 A 的行列式的值 det A 。2. 算法原理2.1 顺序高斯消去法顺序高斯消去法是利用线性方程组初等变换中的一种变换，即用一个不为零的数乘一个方程后加至另一个方程，使方程组变成同解的上三角方程组，然后再 自下而上对上三角

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

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

7、22)La nn)其中m为所进行的行交换次数。这是实际中求行列式值的可靠方法。3.流程图d ?d-bib L,bn 丿开始J从主程序来i+l=i、1ajt, akja,takjj k, k 1,.,nblt,bkbi,tQ1I结束；返回至程床 图3伽选主元屆4. 输出结果- J3事本立电F)精町痕直看蛰賦巴 xCl)=kOO()G x(2)=S. 0000 x：3)=3. 0000detA=-GG- 00005. 结果分析采用计算机运算在计算大数据时有明显的优点，另外也需要考虑到存储。高斯消去法的使用条件是a：）0, k 1,2,L , n，而列选主元法可以保证这一条件。并且可以避免在消元过程

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

9、别化为(A和b(2) ，即M 1 A(1)M b (1)11A(2)b(2)其中m21M1m31Mmn1第 k 次消元时，A(k ) 和 b(k ) 分别化为 A(k 1) 和 b(k 1)M kA(k)M b(k )kA(k1)b(k1)其中mkM1,k Omnk消元过程是对 k 1 n 1 进行的，因此有M n 1 LM 2M 1A(1)A(n)Mn 1 LM 2 M1 (b1) b(n)将上三角形矩阵A(n)记为U,于是有其中L M1 1M 2 1L M1 2 LM11n1M1 21LMn 11Um21m311m32mn1mn2LUm43Mmn3为单位下三角形矩形。这样高斯消去法的实质是

10、将系数矩阵 A 分解为两个三角形矩阵 L 和 U 相乘, 即A LU在上述矩阵描述中遇到了下三角形矩阵运算。主对角线以上元素全为零的方 阵称为下三角形矩阵。下三角形矩阵的乘积仍是下三角形矩阵。若下三角形矩阵 可逆,其逆矩阵仍是下三角形矩阵,而且下三角形矩阵的乘积和逆矩阵很容易求 得。把A分解成一个单位下三角阵和一个上三角阵 U的乘积成为杜里特尔分解。 这种分解是惟一的。2.2 高斯 -约当法高斯消去法有消元和回代两个过程,当对消元过程稍加改变便可以使方程组 化为对角形方程组Dx 的形式,其中矩阵 D 为对角形矩阵,即 a1(11)(2)a22Oa(nnn) nn 的各元素（包括常数项），这个个

11、过程成为归一化，这时方程组的系数阵转化为 单位阵。为减小误差，高斯-约当消去法还常用列选主元技术。当高斯-约当消去法消元的每一步都先用主元去除其所在行的各元素(包括常 数项)时,方程组便可化成11x1x2Mxn(n)b1b2 (n)Mb(n)n这是等号右端即为方程组的解。高斯-约当消去法每一步都用主元去除其所在行3.流程图4.输出结果-CUtpLrt.DCt -记亭萍A的逆矩阵为2.0000 -1.0000 0- 00001. 5000-0.50000- EOOO2. 5000-1.5000 fkSOCX)5. 结果分析采用杜里特尔分解法求解方程组时，由于把对系数矩阵的计算和对右端项的 计算分

12、开，这使计算线性方程组系非常方便。只需进行一次矩形三角分解，然后 再解多个三个方程组，且多解一个方程组仅需要增加大约 n2次乘除法运算。采用高斯约当法仅需要进行消元归一，而不需要回代，为编程实现提供了便 利。6. 心得体会步骤很重要，审题-确定算法-解题步骤-流程图-程序-代入简单值进行验 证。在编程时先在代码输入区打好框架，并且尽量在每一命令后注释，方便检查 错误和日后复习。定义和变量存储很灵活，如我把单位向量直接赋给了A矩阵变量中，还有根据最终的目的直接简化计算。另外赋值前，确定存储空间并且要 定义初值为零。实验四例4-61.题目已知f（X）的观测数据X1234f（X）0-5-63构造插值

13、多项式。2.算法原理首先构造基函数lk（X）Xi，可以证明基函数满足下列条件:其中对于给定n 1个节点,i 0 Xki kl（X）k i1 i kn次拉格朗日插值多项式由下式给出:nL（x）lk（x） ykk 0nX Xi 0 XkXi k由于Ik（X）是一个关于X的n次多项式，所以L（X）为关于X的不高于n次的代数多项式。当X Xi时，L（Xi ） yi，满足插值条件。3.流程图V4.输出结果_ ' Output4.tKt - i3事本交申F)騒帘式心1童看M S期(H)疋二久hOOOOOO处的函数值为：尸-& 3馬00005. 结果分析由于所知的拉格朗日计算机算法只能实现计

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

15、, f (Xi), i 1,2L, m要求在某特定函数类(X)中寻找一个函数(X)作为y f(X)的近似函数，使得二者在Xi上的误差或称残差i(Xi ) f(x ) , i 1,2L, m按某种度量标准为最小，这就是拟合问题。要求残差i按某种度量标准为最小，即要求由残差i构成的残差向量0 1 ,L, mT的某种范数IIII为最小，这本来都是很自然的,为最小。例如，要求maxI 4，或II I即max (X ) f (Xi )可是计算不太方便。所以通常要求:2 3(Xi) f (Xi)i 0m2ii 021或者II2 m2i 0m （X ） f（X ）ii 0为最小。这种要求误差平方和最小的拟合

16、称为曲线拟合的最小二乘法。2.2多变量数据拟合对于给定的一组数据（Xi , yj , i =1, 2,m，寻求做n次多项式nkak X使性能指标1 z 011nJ （a , a 丄，a ）i（yk iX k ）2为最小。k =0,由于性能指标J可以被看做关于ak ,合多项式的构造问题可转化为多元函数的极值问题。令1,n的多元函数，故上述拟J ak从而有正则方程组aoX3XinXn 1XiM2nXia1aMi iMnX yi i对多变量（或称多元）线性模型aoax 1a2x进行了 m次观测%*y2M*yna0a0a1X 11 a2 X21a X a X1 12 2 22a X1 1na X2 2

17、na Xn1 na Xn 2 na Xnn n这个称为回归方程组，写成矩阵形式y=Aa其中1* yyy2* , AM*ymx11x12Mx1mx21x22Mx2mxn1xn2Mxnma0a1Man当 m n 时，要确定一组 ai , i0,1,L, n, 使之精确地满足个方程，这是超定方程组的问题，只能在最小平方误差的基础上确定 i 。定义残差向量 1 , 2 ,L, mT , 则S=y - Aa其中 y y1 , y2 ,L, ym T 代替输出向量。取性能指标使之最小，以此确定出a。由JSTS( y Aa)T ( y Aa) aT A yyT Aa +aT AT Aaa。J = §

18、;T ST=y y利用向量和矩阵的运算公式，有A Aa=A y此即为正则方程组，当 AT A 非奇异时，可求得a= (At A) 1 At y3.流程图4.输出结果9 OLttpLfT5.CKt - 爭本文件旧舷目褂式S1(V)所求参数为: 沪a 2000 bO. 3000 c=0. 50005. 结果分析曲线拟合的最小二乘法是反映所给数据点的总的趋势，并不是严格的通过每 个数据点，这样就避免了大量数据插值时需要高次多项式，同时又去掉了数据所 含的测量误差。6. 心得体会整理思路越来越熟练了，所以执行各个步骤也相对简单了很多。另外对原理 也加深了认识。1#inc_ude=ssio.h=整 nc

19、_ude=3afh.h=#define N 30void 3ain()壬二floaf XLNLGF_LE APUPNfpl Hfopenunpun -Xr-Jr)fp2Hfopen(=oufpun .S-JW-)fscanf(fp 二-%f"cc=fscanf(fpr-cxs 辽裆血fprinff(fp2Jfor(naiAW+)/、卡宜旃并x=+帀x=rx=rQ(2*c*x宇=fprinff(fp2JkHwdux(d)H%.5frrHX=2if(rabs(x=+=-x=DAH0.0§5)break-e-secontinue-fprinff(fp2>n 云w%fHW.3

20、fn?-GX=i2fc-ose(fpDfc_ose(fp2=2> 输入文件： "input1.txt" 18 0.13> 输出文件： “ output1.txt”* 倒数表 * k=0x(0)=0.10000 k=1x(1)=0.06923 k=2x(2)=0.05781 k=3x(3)=0.05564 k=4x(4)=0.05564计算结果：1/18.000000=0.0562. 例 3-41源程序：#include"iostream" #include"cmath" using namespace std;#defin

21、e N 10 void main() int i,j,k,l,n;float bN,aNN,t,d,det=1.0;/* 数据输入 */FILE *fp1,*fp2;fp1=fopen("input2.txt","r");fp2=fopen("output2.txt","w");fscanf(fp1,"%d",&n);for(i=0;i<n;i+) for(j=0;j<n;j+) fscanf(fp1,"%f",&aij);for(i=0;i<

22、;n;i+) fscanf(fp1,"%f",&bi);/* 数据输入 */ /*21*高斯消去 */消元 */列选主元函数 */* /* for(k=0;k<n-1;k+)/ 从第一次消元到第 N-1 次消元d=akk;l=k; for(i=k+1;i<n;i+)/ 找出绝对值最大 的 aik 和 i 行if(fabs(aik)>fabs(d) d=aik;i=l; if(i=n)/ 判断是否奇异，不奇异进行行交换if(d=0) fprintf(fp2," 奇异");/如果所有行的 “首列”都为 0，为奇异 else if(l

23、!=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;i<n;i+)/ 第( k+1 )次消元要得到 N-(k+1) 个乘数aik=aik/akk; for(j=k+1;j<n;j+)/第（i+1）行各列向量对应该行与第（k+1）行各列向量的减法aij=aij-aik*akj;bi=bi-aik*bk;/ 右端常向量/*消元 */* 回代 */bn-1=bn-1

24、/an-1n-1;/ 计算 x(N) 的解for（i=n-2;i>=0;i-）/ 从倒数第二项开始依次回代 N-1 次 t=0;for(j=i+1;j<n;j+)t=t+aij*bj;bi=(bi-t)/aii;/*23数据输出 */* 高斯消去 */* for(i=0;i<n;i+)/ 输出方程组的解fprintf(fp2,"x(%d)=%.4fn",i+1,bi);for(i=0;in;i+) det=det*aii;fprintf(fp2,"detA=%.4fn",det);/ 输出系数矩阵行列式的值fclose(fp1);数据输

25、出 */fclose(fp2);/*2输入文件： “input2.txt ”12 -3 3 -18 3 -115 -15 63输出文件： “output2.txt ”x(1)=1.0000 x(2)=2.0000 x(3)=3.0000 detA=-66.00003. 例 3-101>源程序：#include"iostream" #include"cmath" using namespace std;#define N 30 void main() int i,j,r,k,n;float aNN=0,s;/* 数据输入 */FILE *fp1,*f

26、p2;fp1=fopen("input3.txt","r");fp2=fopen("output3.txt","w");fscanf(fp1,"%d",&n);for(i=0;i<n;i+) for(j=0;j<n;j+) fscanf(fp1,"%f",&aij);for(i=0;i<n;i+)/ 输入单位阵j=i+n;aij=1;/* 数据输入 */*25*LU分解 */for(i=1;i<n;i+)/r=O时："1&qu

27、ot;区间不变化，"2"区间变化。ai0=ai0/a00;for(i=1;i<n;i+)/ 第二行列至末行列进行变化时for(j=i;j<2*n;j+)/ 第"2 (r+1) -1"区间的变化 行 s=0.0;for(k=0;k<i;k+)/ 对应行列元素的乘积进行求和s+=aik*akj;aij=aij-s;for(j=i+1;j<n;j+)/"2 (r+1)" 区间的变化 列 s=0.0;for(k=0;k<i;k+)/ 对应行列元素的乘积进行求和s+=aki*ajk;aji=(aji-s)/aii;

28、/ 现在已将 AI->LUY/* a00=1;/ 第一列其余行已为零*LU分解 */* 高斯约当法解 Ux=Y*/*提 取 UY 减少计算 */for(i=1;i<n;i+)for(j=0;j<i;j+) aij=0;/*提取UY减少计算*/*消元 */for(j=0;j<2*n;j+)/ 首行归一化a0j=a0j/a00;for(i=1;i<n;i+)/(n-1) 次归一消元for(j=i+1;j<2*n;j+)/ 第 i+1 行的各列进行归一化aij=aij/aii;aii=1;for(r=0;r<i;r+)/ 对第一行至第 i-1 行的 i 列进

29、行置零for(j=r+2;j<2*n;j+)/r 行的各列与第 r+1 行的对应列进行减法运算arj=arj-aij*arr+1;/第 r+1 行归一后，乘数即为要置零处的值消元 */arr+1=0;/ 乘数用完之后置零即可/*/*数 据输出 */* 高斯约当法解 Ux=Y*/* fprintf(fp2,"A 的逆矩阵为 nn");for(i=0;i<n;i+) fprintf(fp2,"| ");for(j=0;j<n;j+) fprintf(fp2,"%.4f ",aij+n);fprintf(fp2,"

30、;|n");fclose(fp1);fclose(fp2);/*数据输出 */2> 输入文件： "input3.txt"1 1 -1 1 2 -2 -2 1 13>输出文件： "output3.txt"的逆矩阵为292.0000 -1.00000.00001.5000 -0.50000.50002.5000 -1.50000.50004. 例 4-61源程序：#include"stdio.h" #include"math.h" #define N 50 int n,i,j,k;float xx

31、=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;i<n;i+) fscanf(fp1,"%f,%f",&xi,&yi);/1,0,2,-5,3,-6,4,3 fscanf(fp1,"%f",&xx

32、);/ 要计算的值 for(k=0;k<n;k+)/ 拉 格朗日插值 t=1.0;for(j=0;j<k;j+) t=(xx-xj)*t/(xk-xj);for(j=k+1;j<n;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.txt"1,0 2,-5 3,-6 4,3 2.53>输出文件： "output4.txt&quo

33、t;x=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;/ 定义输入变量float AT3N;/ 定义 A 的转置float X33,Y3;/ 定义中间变量 ATA 和 ATy float s,t,d;FILE *fp1,*fp2;/* 数据输入 */ fp1=fopen("input5.txt","r");fp2=fopen("output5.txt","w");fscanf(fp1,"%d",&n);for(i=0;i<n;i+) fscanf(fp1,"%f,%f,%f",&xi,&yi,&