计算方法实验课程.docx

1、计算方法实验课程计算方法实验课程指导书徐中宇2012年3月目 录第一章 插值法 11拉格朗日插值 12牛顿插值 5第二章 线性方程组的解法 81 高斯消去法 8列主元消去法 123线性方程组的迭代解法 16第三章 方程求根 211 二分法 212 牛顿迭代法 243 埃特金（Aitken）迭代法 26第一章 插值法目的与要求：1. 掌握不同的输入、输出语句，注意节约内存方法。2. 熟悉拉格朗日插值和牛顿插值公式，并体会它们不同的特点。 1拉格朗日插值1方法概要：拉格朗日n次插值多项式 ，其中 可用双重循环来实现；内循环为连乘；外循环为连加。用数组表下标变量，给出。2程序流程图图1-1 拉格朗日

2、插值法程序流程图3程序及例例1已知函数表： xi 0.56160 0.56280 0.56401 0.56521 yi 0.82741 0.82659 0.82577 0.82495用三次拉格朗日插值多项式求x=0.5635时的函数值。程序为#include stdafx.h#include const int n=3;const double xx=0.5635;void main() double xn+1=0.56160,0.56280,0.56401,0.56521; double yn+1=0.82741,0.82659,0.82577,0.82495; double L=0,Li=

3、1; for (int i=0;i=n;i+) Li=1; for (int j=0;j=n;j+) if (j!=i) Li=Li*(xx-xj)/(xi-xj); L+=Li*yi; coutThe result is:Lendl;例2 已知函数表，用拉格朗日插值多项式求0.5，0.7，0.85三点处的函数值。程序的结构：在程序的说明部分将已知数据分别存入一维数组x和y。在主函数main中提示人工键入插值点的数目m，继而在m次循环计算中分别提示键入插值点，计算插值，并输出结果。程序中的主要常量，变量及函数：n为插值基点的最大下标。xn+1,yn+1分别存放插值基点及其函数值。M为插值点个数

4、，X存放插值点的值。程序清单：#include stdafx.h#include #define n 4 /*插值基点的最大下标*/void main() double xn+1=0.4,0.55,0.8,0.9,1; /*插值基点的值*/ double yn+1=0.41075,0.57815,0.88811,1.02652,1.17520; int m,k,i,j; float X,L,P; printf(n Please enter m=); scanf(%d,&m); /*键入插值点的个数*/ for(k=1;k=m;k+) printf(n Please enter X%d=,k);

5、 scanf(%f,&X); P=0.0; for (i=0;i=n;i+) L=1.0; for (j=0;j=n;j+) if (j!=i) L=L*(X-xj)/(xi-xj); P=P+yi*L; printf(P(%f)=%fn,X,P); 计算结果：（当提示“Please enter m=”时，键入插值点个数3，再分别根据提示X1=或X2=，X3=依次键入插值点0.5,0.7,0.85），输出的结果分别为P（0.500000）=0.521090 P(0.700000)=0.758589 P(0.850000)=0.9561192牛顿插值1方法概要：牛顿n次插值多项式 Pn(x)=P

6、n-1(x)+an(x-x0)(x-x1)(x-xn-1) 在计算机上计算均差表时，可用一维数祖x存放基点值，用二维数组F存放和计算由个阶均差构成的下三角矩阵。计算公式如下：Fi0=f(xi) ,i=0,1,n Fij=(Fi,j-1-Fi-1,j-1)/(xi-xi-j) ,j=1,2,n i=j,j+1,n2. 程序及例例3. 已知函数表同例2，试根据上面公式的算法计算均差表存入矩阵F，并用牛顿均差插值多项式求解上题。注意到，那么 计算的流程为：（1）。 （2）对，计算。程序的结构：在程序开始时将初始数据分别存入一维数组x,y。主函数main首先调用自定义函数junca计算均差表，再提示人

7、工键入插值点个数m，继而在m次循环中分别提示键入插值点X，并计算相应的插值，输出值结果。程序中的主要常量，变量及函数：n,x,y,m,X均同上例。Junca(double x,double y,double Fn+1)是自定义函数，计算均差表并打印输出，结果存于F中，P用于计算牛顿均差插值多项式的值。#include stdafx.h#include const int n=4; /*函数junca实现均差表的计算，将均差表存入数组F，并在屏幕上打印输出*/void junca (float x,float y,float Fn+1) int i,j; for (i=0;i=n;i+) Fi0

8、=yi; for (j=1;j=n;j+) for (i=j;i=n;i+) Fij=(Fij-1-Fi-1j-1)/(xi-xi-j); printf (n %12s%12sn,Xi,F(Xi); for (j=1;j=38;j+) printf(-); printf (n); for (i=0;i=n;i+) printf(%12f,xi); for (j=0;j=i;j+) printf(%12f,Fij); printf(n); for (j=1;j=38;j+) printf(-); printf(n);void main() float xn+1=0.4,0.55,0.8,0.9,

9、1; /*插值基点的值*/ float yn+1=0.41075,0.57815,0.88811,1.02652,1.17520; int i,m,k; float p,t,Fn+1n+1; float X; junca (x,y,F); /*调用函数junca计算均差表*/ printf(n Please enter m=); scanf(%d,&m); /* 键入插值点的个数 */ for (i=1;i=0;k-) p=p*(X-xk)+Fkk; /*计算牛顿均差值多项式*/ printf(P(%f)=%fn,X,p); 计算结果：XiF(Xi)1阶2阶3阶4阶0.4000000.4107

10、500.5500000.5781501.1160000.8000000.8881101.2398400.3096000.9000001.0265201.3841000.4121710.2051431.0000001.1752001.4868000.5135000.2251750.033386(当提示Please enter m=时，键入插值点个数3，再分别根据提示X1=或X2，X3=依次键入插值点0.5,0.7,0.85)输出的结果为P0.500000=0.521090 , P0.700000=0.758589 , P0.850000=0.956119练习题1 已知函数表：xi 0.4 0.5

11、5 0.65 0.8 0.9 yi 0.41075 0.57815 0.69675 0.88811 1.02652分别用四次拉格朗日插值和牛顿插值，求x=0.596时的函数值。第二章 线性方程组的解法本章的目的与要求： 1 熟悉求解线性方程组的有关理论与方法。 2 会编制列主元消去法、LU分解法与两种迭代法的计算程序，上机计算求得结果， 3 通过上机实习对计算数据的分析，进一步了解各种方法的功能，优缺点与适用范围，体会如何针对不同问题适当选择方法。1 高斯消去法 1. 方法概要 解线性方程组 (1)为统一计算公式与简化计算程序，把方程组的增广矩阵记作A(1)=(A,b) (2)解题分为消元与回

12、代两个过程：消元过程；第1步，对第1列作消元，将消为零。第2步，对第2列作消元，将消为零。设0(k=1,2,3,4.,n-1)。一般地,对于第k步消元,采用如下公式: (3) (4) (j=k,k+1,n+1)直到作完-1步消元，相应的方程组化为与原来方程组等价的方程组： （6） 回代过程：回代过程是解方程组，采用如下公式 （7） 2程序流程图(见图21) 流程图说明(1)首先输入方程组的阶数n及方程组的增广矩阵。为一个(+1)的二维数组：所有运算的中间结果与最后结果全部放中相应的位置。(2)本程序为了教学的需要，先屏幕输出原始增广矩阵，然后逐次输出每次消元后的矩阵以及每次所用的消元因子，最后

13、输出解向量。我们可以从中看出消元过程的具体情况。在实用中无此必要，可以省略。(3)这是顺序取主元的消去法。为了检查是否为零，设置了一个检查口，一旦遇到，便使其输出一个记号，以示顺序取主元的方法失败. 3程序及例 例1.求解方程组 2.5x+2.3x-5.1x3= 3.7 5.3x+9.6x+1.5x= 3.8 8.lx+1.7x2-4.3x=5.5(答：x=(0.406705, 0.237818,-0.419325)T)#include stdafx.h#include #include #define n 3 /*方程组的阶数*/double aann+1=2.5,2.3,-5.1,3.7,

14、5.3,9.6,1.5,3.8,8.1,1.7,-4.3,5.5; /*增广矩阵的初始数据*/int gauss (double a n+2,double x ) int i,j,k; float c; for (k=1;k=n-1;k+) /*消元过程*/ for (i=k+1;i=n;i+) if (fabs(akk)1e-12) printf(n det=0.fail!n); return(0); c=aik/akk; for (j=k+1;j=1;k-) /*回代过程*/ xk=akn+1; for(j=k+1;j=n;j+) xk=xk-akj*xj; xk=xk/akk; retu

15、rn(1);void main() int i,j,det; double an+1n+2,xn+1; for(i=1;i=n;i+) for(j=1;j=n+1;j+) /*用a11ann+1存放增广矩阵*/ aij=aai-1j-1; det=gauss(a,x); /*调用函数gauss求解方程组，并获取返回标志值*/ if (det!=0) for(i=1;i=n;i+) printf(nx%d=%f,i,xi); printf(n-n); 是 图21 Gauss消去法程序流程图列主元消去法1方法概要作第K步消元时，从A（K）的子矩阵（见1中式（3）的第一列中选取, 以作为主元素，第行

16、称为主行。将主行与第k行交换，然后再作这一步的消元。2计算框图 图 22 列主元消去法框图3程序及例例2.用列主元消去法求解 程序的结构：主函数main调用函数gaussl进行列主元高斯消去法计算，根据返回的函数值为1或0而输出计算结果或失败信息。程序中的主要常量、变量及函数：n为常量，指定方程组的阶数。数组aann+1存放方程组增广矩阵的原始数据并保留不变。数组an+1n+2、xn+1分别用于计算增广和方程组的解。int gauss1(a,x)是自定义函数，用列主元高斯消去法求解n阶线性方程，当系数矩阵奇异使求解失败时，输出信息“det=0.fail!”并返回函数值0；否则求解成功，返回函数

17、值1。Int det在主函数里接受gaussl的返回值，当det不为0时主函数打印出方程组的解。#include stdafx.h#include #include #define n 3 /*方程组的阶数*/static double aann+1=1,2,-1,3,1,-1,5,0,4,1,-2,2; /*增广矩阵的初始数据*/int gauss1 (double a n+2,double x ) int i,j,k,r; double c; for (k=1;k=n-1;k+) /*消元过程*/ r=k; for (i=k;ifabs(ark) r=i; if (fabs(ark)1e-

18、12) printf(n det=0. fail! n); return(0); if (r!=k) for (j=k;j=n+1;j+) /*交换k、r两行*/ c=akj; akj=arj; arj=c; for (i=k+1;i=n;i+) /*进行消元计算*/ c=aik/akk; for (j=k+1;j=n+1;j+) aij=aij-c*akj; if (fabs(ann)=1;k-) /*回代过程*/ xk=akn+1; for(j=k+1;j=n;j+) xk=xk-akj*xj; xk=xk/akk; return(1);void main() int i,j,det; d

19、ouble an+1n+2,xn+1; for(i=1;i=n;i+) for(j=1;j=n+1;j+) /*用a11ann+1存放增广矩阵*/ aij=aai-1j-1; det=gauss1(a,x); /*调用函数gauss1求解方程组，并获取返回标志值*/ if(det!=0) for(i=1;i=n;i+) printf(nx%d=%f,i,xi); printf(n-n); 计算结果： x1=0.250000 x2=1.500000 x3=0.2500003线性方程组的迭代解法1. 雅可比迭代法例3. 根据雅可比迭代法编制程序求解方程组。其中的计算流程为：（1）;（2）对做： 若

20、则。程序中的主要常量，变量：n，eps是常量，表示方程组的阶数及指定的精度。aann,bbn分别存放方程组的系数矩阵及常数项。int NO是最大迭代次数，当迭代达到NO次以上时程序输出失败信息，并结束。xn+1,yn+1分别用于存放相继两次的迭代向量。用变量max计算。程序清单:/ Defines the entry point for the console application./#include stdafx.h#include #include #define n 3 /*方程组的阶数*/#define eps 0.5e-4 /*给定精度要求*/static double aann=

21、10,-1,-2,-1,10,-2,-1,-1,5;static double bbn=7.2,8.3,4.2;void main() int k=0,NO,i,j; double d,s,max; double an+1n+1,bn+1,xn+1,yn+1; for (i=1;i=n;i+) for (j=1;j=n;j+) aij=aai-1j-1; bi=bbi-1; printf(n Please enter No:); scanf(%d,&NO); /*键入最大迭代次数NO*/ for (i=1;i=n;i+) yi=0; do k+; for (i=1;i=NO) printf(n

22、 Fail !); break; /*当k=NO时输出失败信息，结束迭代*/ max=0.0; for (i=1;i=n;i+) s=0.0; for (j=1;j=n;j+) if(j!=i) s=s+aij*xj; yi=(bi-s)/aii; d=fabs(yi-xi); if (max=eps); if (maxeps) /*当maxeps时结束迭代*/ printf(k=%dn,k); for (i=1;i=n;i+) printf(x%d=%fn,i,xi); 运行结果：产生输入提示Please enter NO：后，键入适当大的正整数，可行到结果数据为 k=11 x1=1.099

23、979 x2=1.199979 x3=1.2999752. Gauss-Seidel迭代法例4. 根据高斯-赛德尔迭代法编制程序求解上题的方程组。程序中的主要常量及变量：n,eps,aa,bb,a,b,x，NO均同上例。每次迭代计算时，用变量s临时存放xi的旧值。用d计算xi与其旧值之差最大绝对值。程序清单:/ Defines the entry point for the console application.#include stdafx.h#include #include #define n 3#define eps 0.5e-4static double aann=10,-1,-2

24、,-1,10,-2,-1,-1,5;static double bbn=7.2,8.3,4.2;void main() int k=0,NO,i,j; double d,sum,s; double an+1n+1,bn+1,xn+1; for (i=1;i=n;i+) for (j=1;j=n;j+) aij=aai-1j-1; bi=bbi-1; printf(n Please enter NO:); scanf(%d,&NO); for (i=1;i=NO) printf(n Fail !); break; for (i=1;i=n;i+) sum=0.0; s=xi; /*s临时存放xi的旧值*/ for (j=1;j=n;j+) if(j!=i) sum=sum+aij*xj; xi=(bi-sum)/aii; /*计算xi的新值*