数值计算方法

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) 。将非线性方程 f x( ) = 0的根 x*代入 f x( *) = 0 ，即f x( k ) + f (xk )(x* xk ) 0*xk f (xk ) 解出x f (xk )将右端取

2、为 xk+1 ，则 xk+1 是比 xk 更接近于 x* 的近似值，即f (xk )xk+1 xk f (xk ) 这就是牛顿迭代公式，相应的迭代函数是f (x)(x) = x f (x)2.2 牛顿迭代法的应用11计算 是求cx =1 0的解，解出 x，即得到 。取 cc 有牛顿迭代公式cxk 11 xk+1 = xk = cc 这样就失去了迭代的意义，达不到迭代的效果。1f (x) = cx1， f (x) = c，故重新构造方程： cx2 x = 0 ， 也是该式的解。故取 f (x) = cx2 x ，cf (x) = 2cx 1，则有牛顿迭代公式xk+1 = xk cxk2 xk =

3、 cxk2,k = 0,1,.2cxk 12ck 111的值在 之间，取初值 x0 = 0.1。20103.流程图0,Nx读入1k()0?0xf=1x输出011kkxx+()()0100fxxxfx10?xx=4.输出结果5.结果分析当k= 3时，得 5 位有效数字 0.05 564。此时， x3 x4 = 0.00 000 x* 时必收敛，但是当 x0 0) 时迭代结果发散，较小尚不确定。6.心得体会起初对题目的理解并不是很透彻，另外对构建牛顿迭代公式理论依据不是特别充分，比如说为什么在原有直接得到的式子两边各乘一个 x，只是试出来的。在编程方面不够成熟。当然也加深了对牛顿迭代法的理解和应用

4、的具体实现。实验二 例 3-41.题目用列主元消去法求解方程组 12x1 3x2 + 3x3 =1518x1 3x2 x3 = 15 x1 + x2 + x3 = 6 并求出系数矩阵 A的行列式的值det A。2.算法原理2.1 顺序高斯消去法顺序高斯消去法是利用线性方程组初等变换中的一种变换，即用一个不为零的数乘一个方程后加至另一个方程，使方程组变成同解的上三角方程组，然后再自下而上对上三角方程组求解。这样，顺序高斯消去法可分成“消去”和“回代” 两个过程。在用顺序高斯消去法时，在消元之前检查方程组的系数矩阵的顺序主子式，当阶数较高时是很难做到的。若线性方程组的系数具有某种性质时，如常遇到的

5、对角占优方程组，自然能够用高斯消去法求解。2.2 列选主元消去法线性方程组只要系数矩阵非奇异，就存在惟一解，但是按顺序消元过程中可能出现主元素akk( )k = 0，这时尽管系数矩阵非奇异，消元过程无法再进行，或者即使akk( )k 0，但如果其绝对值很小，用它作除数也会导致其他元素的数量级急剧增大和使舍入误差扩大，将严重影响计算的精度。为避免在校园过程确定乘数时的所用除数是零或绝对值小的数，即零主元或小主元，在每一次消元之前，要增加一个选主元的过程，将绝对值大的元素交换到主对角线的位置上来。列选主元是当高斯消元到第 k 步时，从 k 列的akk 以下（包括akk ）的各元素中选出绝对值 大的

6、，然后通过行交换将其交换到akk 的位置上。交换系数矩阵中的两行（包括常数项），只相当于两个方程的位置交换了，因此，列选主元不影响求解的结果。列选主元消去法常用来求行列式。设有矩阵a11 L a1n A = MM an1 L ann 用列主元消去法将其化为上三角形矩阵，对角线上元素为a11(1) ,a22(2) ,L,a33(3),于是行列式det A = (1)ma a11(1) 22(2) Lann( )n其中 m 为所进行的行交换次数。这是实际中求行列式值的可靠方法。3.流程图4.输出结果5.结果分析采用计算机运算在计算大数据时有明显的优点，另外也需要考虑到存储。高斯消去法的使用条件是a

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

8、数矩阵 A进行除不交换两行位置得初等行变换相当于用初等矩阵M1 左乘 A，在对方程组第一次消元后，A(1) 和b(1) 分别化为 A(2) 和b(2) ，即M1 A(1) = A(2)M1b(1) = b(2) 1m211其中M1 = m311 MOmn11第 k 次消元时， A( )k 和b( )k 分别化为 A(k+1) 和b(k+1) ，即M Ak ( )k = A(k+1) M b( )k = b(k+1)k1其中 M1 = O1mk+1,k MmnkOO1消元过程是对k=1 n1进行的，因此有Mn1LM M A21(1) = A( )nM LM M b2(1) = b( )n n 1

9、1 将上三角形矩阵 A( )n 记为 U，于是有A = M M1121LMn11U = LU其中 1m211mm1L = M M1121LMn11 = 3132 MMm43 O MMMOmn1mn2mn3 L 1为单位下三角形矩形。这样高斯消去法的实质是将系数矩阵 A 分解为两个三角形矩阵 L 和 U 相乘，即A=LU在上述矩阵描述中遇到了下三角形矩阵运算。主对角线以上元素全为零的方阵称为下三角形矩阵。下三角形矩阵的乘积仍是下三角形矩阵。若下三角形矩阵可逆，其逆矩阵仍是下三角形矩阵，而且下三角形矩阵的乘积和逆矩阵很容易求得。把 A 分解成一个单位下三角阵和一个上三角阵 U 的乘积成为杜里特尔分

10、解。这种分解是惟一的。2.2 高斯-约当法高斯消去法有消元和回代两个过程，当对消元过程稍加改变便可以使方程组化为对角形方程组Dx b= 的形式，其中矩阵D为对角形矩阵，即a11(1)(2)D = a22Oann( )n 当高斯-约当消去法消元的每一步都先用主元去除其所在行的各元素（包括常数项）时，方程组便可化成1 x1 b1( )n 1 O xM2 = bM2( )n 1 xn ( )n bn这是等号右端即为方程组的解。高斯-约当消去法每一步都用主元去除其所在行的各元素（包括常数项），这个个过程成为归一化，这时方程组的系数阵转化为单位阵。为减小误差，高斯-约当消去法还常用列选主元技术。4.输出

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

12、验四 例 4-61.题目已知 f (x) 的观测数据x1234f( )x0-5-63构造插值多项式。2.算法原理首先构造基函数l xk ( )=i=n0 xxk xxii ，可以证明基函数满足下列条件：i k0i kl xk ( i ) = ， 1 i = k对于给定n+1个节点，n次拉格朗日插值多项式由下式给出：nL x() = l x yk ()kk=0其中l xk ( )=i=n0 xxk xxiii k由于l xk ( ) 是一个关于 x的n次多项式，所以L x( ) 为关于 x的不高于n次的代数多项式。当 x = xi 时，L x( i ) = yi ，满足插值条件。3.流程图1t0

13、0yk,n0,1,iii=L输入（x,y）0,k1,k1,njkjxxttxxj+=yktyy+=4.输出结果5.结果分析由于所知的拉格朗日计算机算法只能实现计算某一特定值的近似函数值，而不知如何导出表达式，故例求 x=2.5 处的函数值以说明表达式以得出，只是在计算机程序中。并且也能达到拉格朗日插值法使用的目的。6.心得体会编程不够细心，程序没问题，却因为不知道是输入文件错了而检查了好长时间。但同时也加深了对拉格朗日基函数性质的认识和理解。实验五 习题 5-21.题目给出平面函数 z x y( ,) = ax +by +c的数据i12345xi0.10.20.40.60.9yi0.20.30

14、.50.70.8zi0.580.630.730.830.92按 小二乘原理确定a b c， 。2.算法原理2.1 小二乘原理设已知某物理过程 y = f x( )的一组观测数据(xi , f x( i ),i =1,2,L,m要求在某特定函数类( )x 中寻找一个函数( )x 作为 y = f x( )的近似函数，使得二者在 xi 上的误差或称残差i =(xi ) f x( i ) , i =1,2,L,m按某种度量标准为 小，这就是拟合问题。1)ii要求残差i 按某种度量标准为 小，即要求由残差i 构成的残差向量 =0, 1,L, m T 的某种范数 为 小。例如，要求 ，或 即mm 1 =

15、 i = (xi ) f x( i=0i=0 = max i = max ( )xi f x( )ii为 小，这本来都是很自然的，可是计算不太方便。所以通常要求：11 2 =m i2 2 =m ( )xi f x( )i 2 2 i=0 i=0或者mm 22= i2 =(xi ) f x( i ) 2i=0i=0为 小。这种要求误差平方和 小的拟合称为曲线拟合的 小二乘法。2.2多变量数据拟合对于给定的一组数据 (xi , yi ) ，i=1，2，m ，寻求做n次多项式ny =a xkkk=0使性能指标J a a( 0,1,L,an ) = (yi a xkik )2 为 小。i=1k=0由于

16、性能指标J可以被看做关于ak ，k=0，1，n的多元函数，故上述拟合多项式的构造问题可转化为多元函数的极值问题。令J = 0从而有正则方程组 m xixi23 xixi2xi MMMLLLxin a0 yi n+1 a1 xi yi xi M = MM ak xinxin+1xin+2 L xi2n anxin yi 对多变量（或称多元）线性模型y* = a0 +a x1 1 +a x22 +L+a xnn进行了m次观测y1* = a0 +a x1 11 +a x221 +L+a xn1ny2* = a0 +a x1 12 +a x222 +L+a xn2n Myn* = a0 +a x1 1

17、n +a x22n +L+a xnnn这个称为回归方程组，写成矩阵形式y=Ay1* 1x11 * 其中y= 2 , A=1x12yM MM * 1x1mym x21 L xn1 a0 x22 L xn2 , =a1 。M M M M x2m L xnm an 当m n 时，要确定一组a ii , = 0,1,L,n, 使之精确地满足 m 个方程，这是超定方程组的问题，只能在 小平方误差的基础上确定i 。定义残差向量= 1,2,L,m T , 则= y-A其中 y = y y1,2,L, ym T 代替输出向量。取性能指标J = T使之 小，以此确定出。由J = T = (y A) (T y A

18、) =y yT TAy yTA+TA AT 利用向量和矩阵的运算公式，有A AT =AT y此即为正则方程组，当A AT 非奇异时，可求得 = (A AT)1AT y3.流程图nizyxiii,2,1,L=输入n5niATyATxATiiiii,.,2,1;1321=k1kniikikjnijikiYzajXaa=113,2,13?=kkk+1=输出参数4.输出结果5.结果分析曲线拟合的 小二乘法是反映所给数据点的总的趋势，并不是严格的通过每个数据点，这样就避免了大量数据插值时需要高次多项式，同时又去掉了数据所含的测量误差。6.心得体会整理思路越来越熟练了，所以执行各个步骤也相对简单了很多。另

19、外对原理也加深了认识。附录：1.造倒数表1源程序：#includestdio.h#includemath.h #define N 30 void 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)=%.

20、5fn,i,i,xi); if(fabs(xi+1-xi)输入文件：input1.txt180.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.=0.0562.例 3-41源程序：#includeiostream #includecmath using namespace std; #define N 10 void main()int i,j,k,l,n; float bN,aNN,t,d,det=1.0;/*数据

21、输入*/FILE *fp1,*fp2; 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)/判断是否奇异，不奇异进行行

22、交换 if(d=0)fprintf(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=(b

23、i-t)/aii;/* *高斯消去*/*数据输出*/ for(i=0;in;i+)/输出方程组的解fprintf(fp2,x(%d)=%.4fn,i+1,bi); for(i=0;i输入文件：“input2.txt”312 -3 3-18 3 -11 1 115 -15 63输出文件：“output2.txt”x(1)=1.0000 x(2)=2.0000 x(3)=3.0000 detA=-66.00003. 例 3-101源程序：#includeiostream #includecmath using namespace std; #define N 30 void main()int i

24、,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=1;/*数据输入*/*LU 分解*/ for(i=1;in;i+)/r=0 时：1区间不变化，2区间变化。ai0=ai0/a00;for(i=1;in;i+)/第二行列至末行列进行变化时 for(j=i

25、;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*/*提取 UY 减少计算*/for(i=1;in;i+) for(j=0;ji;j+) aij=0;/*提取 UY 减少计算*/*消元*/ for(j=0;j2*n;j+)/首行归一化a0j=a0j/a00;a00=1;/第一列其余行已为零for(i=1;in;i+)/(n-1)次归一消元 fo

26、r(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;/乘数用完之后置零即可/*消元*/*高斯约当法解 Ux=Y*/*数据输出*/ fprintf(fp2,A 的逆矩阵为nn);for(i=0;in;i+) fprintf(fp2,|); for(j=0;j输入文件：input3.txt31

27、1 -11 2 -2-2 1 13输出文件：output3.txtA 的逆矩阵为|2.0000-1.00000.0000|1.5000-0.50000.5000|2.5000-1.50000.5000|4.例 4-61源程序：#includestdio.h#includemath.h#define N 50int n,i,j,k; float xx=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;j输入文件：input4.txt41,02,53,64,32.53输出文件：output4.txtx=2. 处的函数值为:y=-6.5.习题 5-21源程序：#include stdio.h#include math.h #define N 30 void