数值分析综述报告

淮阴工学院数值分析考试基于 Matlab 的方法综合应用报告班 级： 金融1121名：学 号：成 绩：数理学院2014年6月7日数值分析课程综述报告前言： 数值分析也称计算方法，它与计算工具的发展密切相关。数值分析是一门为科学计算提 供必需的理论基础和有效、实用方法的数学课程，它的任务是研究求解各类数学问题的数值 方法和有关的理论。正文：第一章 近似计算与误差分析模型误差；观测误差；截断误差；舍入误差。1、产生误差的原因2、四则运算的误差 加减法运算8 C* y* )= 8 C* )+6 (y*)乘法运算xy x * y * | = | lx|y y *n 6xy - xy *y (x* y * ) = |x+ xy * x * y * |x - x* I)6 (y * )+ | y *(6 (x* )xx *xy * x * yyy *yy *xy * x * y 除法运算：*yy(x - x * ) y * + ( y * - y ) x *yy *6 ( x * ) I y * I(y3、科学表示法、有效数字、近似值的精度 任何一个实数都可以表示成如下的形式：其中：是正整数，是整数，如果是数的近似值并且则称该近似值具有位有效数字(significant digit)。此时，该近似值的相对误差为另一方面，若已知那么，6r (* ) 占101n1x x*| |x*|6 C*)r10m1 x a .a a “101-n2 (a +1)1 丄 10m n2即：x*至少有n位有效数字。例：兀=3.141592653589793 取其近似值如下： x*=3.14 x*=3.14159 x*=3.1415 x*=3.141x* = 10 x 0314b-x*| = o.o。16 o.。05 =1 x10-2 = 2x1013x* = 10 x 0.314159|兀x*| = 0.0000026. 0.000005 =丄 x10-5 =丄 x101-62 2x* =10x0.31415n x*| = 0.000092. 0.0001 丄 x 10-3 =丄x 101-42 2x* = 10x 0.3141兀x*=0-00059- 0-001 2x10 -2=2x101-3第二章 线性方程组 在科学计算中，问题的本身就是求解线性方程组，许多问题的求解需要最后也归结为线 性方程组的求解，所以线性方程组的求解是科学计算中最常见的问题。对于线性方程组的求解一般有两种方法：（1） 直接法：高斯消去法；（2） 间接法：各种 迭代法。（1） 高斯消去法： 求解思路：先消元，即按一定的规律逐步消去未知量，将方程组Ax = B化为等价的上（或下）三角形方程组；然后进行回代，即由上三角形方程组逐个求出； 高斯（列、全）主元素消去法及在消元的每一步选取（列）主元素一一列中绝对值最 大的元取做主元素，计算步骤：消元过程：按列选主元、行交换、消元计算；回代过程； 高斯列主元素消去法的MATLAB实现：。第三章解线性方程组的迭代法通常逆矩阵不易求得，特别是对于大型的线性方程组，需要用迭代法求解。 用迭代法求解线性方程组，要把线性方程组写成等价的形式，右边写为迭代格式，如Ax = -b n x = Mx + b n x = Ax + x + b = （A + E） x + b = Mx + bn kx = Ax + kx + b = （A + k E ） x + bnx = 1 （A + k E ）x + kn k2、关于迭代法收敛性的两个重要结论： 充分必要条件是：矩阵的谱半径Z、；p（M ） 1 充分条件是：矩阵旳的某个算子范数“ “。Ml 1 13、线性方程组的迭代法主要有Jacobian迭代法，Gauss-Seidel迭代法。Jacobian迭代法：Ax = b/、qn( D L U ) x = b、A = D L UnDx = ( L + U ) x + bnx (k + 1) = D i ( L + U ) x (k)+ D i b M = D - i ( L + U )nq、f = D -1 bGauss-Seidel 迭代法：Ax = b/、qn (D L ) x = Ux + bA = D L Ux(k+1)=(D L)1Ux(k)+(D L)1bnM =(DL)1Uf = (D L )-1 b(3.7)Jacobian迭代法与G-S迭代法比较:x (k+1)、1X (k+1)2.一 000X (k+1)、1X (k+1)2.=D-1a2100X (k+1) 匕丿an1a0nn1j X (k+1)?n-0a-12. a.Inf X (k 八100 X (k )J一an1n2_ 00 0 _X (k ) 丿+ D-ib+D-1n式（3.7）和（3.8）表明：Gauss-Seidel迭代法在计算第k +1次迭代的第i个分量x（k+i）时，i及时地利用了在此步迭代中得到的新的迭代值：皿+1），j = ,2,-，由于第k +1步的迭 代值通常比第k步的迭代值更接近方程组的精确解，所以，在Jacobian迭代法和GS迭代法 都收敛的情况下，Gauss-Seidel迭代法的收敛速度比Jacobian迭代法的收敛速度高。例题：用MATLAB函数normrdn生成5阶矩阵M和向量b分别构造线性方程组Ax = b的 Jacobi迭代格式和G-S迭代格式，并判断收敛性。Jacobian迭代法和GS迭代法程序如下：clc;clear all;%1、生成M和bM=normrnd(1,2,5) b=normrnd(1,2,5,1) %Jacobian迭代法M1=D(L+U) f1=Db rho=max(abs(eig(M1);R=1e-08;%设定的一个收敛标准switch sign(1-rho)case -1disp(the Jocobian method is not applicable) otherwisex(:,1)=normrnd(0,9,5,1);k=1while k=R k=k+1;elseX=x(:,k+1);disp(Jacobian迭代法迭代次数为：) IterN=k %Jacobia n迭代法迭代次数 breakendendend%Causs-Seide 1 迭代法M2=(D-L)U f2=(D-L)b rho=max(abs(eig(M2);R=1e-08;switch sign(1-rho)case -1disp(the auss-seide1 method is not app1icab1e)otherwise x(:,1)=normrnd(0,9,5,1); k=1 while k=R k=k+1; else X=x(:,k+1) disp(Causs-Seidel迭代法迭代次数为：) IterN=k breakendendend第四章 非线性代数方程(组)的数值解法：一、二分法：首先要确定适当的包含根的区间，这可以依据闭区间上连续函数的介值定理来确定，例如该方程：f (x) =x 3 一 x 2 一 0.8 x + 0.75 = 00对于该方程f (-1)=-2 + 0.8 + 0.75 0所以该方程至少有一个实根位于区间，图像表明该区间中只含有一个实根；用x*表示方程f(x)= 0在区间a,b上的精确解，对于给定的精度要求* 0,取区间切 的中点，并按下式进行判断：a + bx =2 ( fix 丿=0 n x = x *f (x) f(a) 0 nx * e a ,x f (x) f(b) 0 nx * e x ,b 以 f (x )f (a) 0 为例, 如果没有达到精度要求，令X = b ,并重复(1)的迭代过程；ba2如果 j，那么，必有|x-x* , Vx轟x2 =Frtic;f(:,k)=F(X(1,k),X(2,k),X(3,k);J=subs(DFx,X(1,k);J=subs(Jy,X(2,k);J=subs(J,z,X(3,k);X(:,k+1)=X(:,k)_Jf(:,k);t(k)=toc;if norm(X(:,k+1)_X(:,k),2)Erbreakendk=k+1;enddisp( Newton迭代法结果为：);disp(X(:,end);运行结果：Newton 迭代法结果为：3.42910.46580.6535第五章插值一、插值：插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况, 估算出函数在其他点处的近似值。二、常用差值法：拉格朗日(Lagrange)插值法、牛顿(Newton)插值法1、拉格朗日(Lagrange)插值法：拉格朗日插值多项式L(x)=工 L (x)ykkk 二 0简单的证明因为拉格朗日插值多项式的基函数有如下的性质：k = 0,1,2,nl (x)=n kxxkx x , j H kj所以拉格朗日插值多项式k 0,1,2,nL(x )=X L (x )y L (x )y ykj k j k k k kj=0满足插值的条件。插值多项式L(x)= Lk (x)yk拉格朗日插值法的不足在实际问题中，观测的数据可能会不断增加，如果用拉格朗日插值公式构造插值多项 式，那么，每当增加数据就要重新计算多项式的系数，由此增加许多不必要的计算工作量。2、(三次)样条(Spline)插值插值条件要求要求所求的插值多项式S (x)(三次样条函数) 在每个区间k一 1x，k = 0，1，2,n，是次数不超过三次的多项式; S(x )= y , k = 0,1,2 * *, n. S(x)在区间上具有二阶连续导数。例题：在夏季，较大湖泊的水体按深度被跃变层分为上部的变温层和下部的均温层。水体的Zt (x )dx分层化对环境工程中污染问题的研究具有重要的意义，例如，有机物的分解将导致被跃变层 隔离的底部水体中氧急剧减少。按温度随水深的变化曲线T - T(x)，跃变层位于水深处：深度(m )： x02.34.99.113.718.322.927.2温度(0C): T22.822.822.820.613.911.711.111.1现有美国普拉特湖0max0 x bottom(Platte Lake)的一组数据:1.试利用Lagrange差值方法求温度随水深(近似)变化函数卩=Tl(x)表达式;2.试利用三次样条差值方法(应用Mat lab函数csape)求温度随水深(近似)变化函数T T (x)表达式.s3.画出插值函数T = T (x)和T T (x)曲线，并与原始插值数据图像作比较 ls程序代码：函数文件程序：function Ln = my_Fun(x, XI, YI)if isa(x, sym) = 1;n = length(XI) - 1;Ln = 0;Pn = sym(ones(n + 1, 1);for k = 1 : n + 1XI(i);XI(i);for i = 1 : n + 1 if i = kPn(k) = Pn(k)*(x - XI(i)/(XI(k) elsePn(k) = Pn(k);endendLn = Ln + Pn(k)*YI(k);endelsen = length(XI) - 1;L = ones(n + 1, length(x);Ln = zeros(size(x);for k = 1 : n + 1for i = 1 : n + 1if i = kL(k,:) = L(k,:).*(x - XI(i)./(XI(k) elseL(k,:) = L(k,:);endendLn = Ln + L(k, :).*YI(k);endEnd主文件程序：clc;clear all;close all;%问题1用Lagrange差值方法求温度随水深(近似)变化图像t = sym(t, real);x=0 2.3 4.9 9.1 13.7 18.3 22.9 27.2;T=22.8 22.8 22.8 20.6 13.9 11.7 11.1 11.1;X=linspace(0,27.2,275);Ln=my_Fun(X,x,T);figure(1)plot(x,T,r*,markersize,15)xlabel(深度x);ylabel（温度T）;title（ 原始的散点图）;pause（3）holdonplot（X,Ln,b-,linewidth,3）;xlabel（深度x）;ylabel（温度T）;title（ Lagrange差值图）;%求出温度随水深（近似）变化函数表达式I_poly = my_Fun（t, x, T）;I_poly = simple（I_poly）I_poly = sym2poly（I_poly）%2.用三次样条差值方法求温度随水深（近似）变化函数 表达式figure（2）holdon y=interp1（x,T,X,spline）;p l o t （ X , y , b- , l i n ewidth , 3 ） ; %用三次样条插值方法画出图像 xlabel（深度x）;ylabel（温度T）;title（ 三次样条法差值图）;PP=csape（x,T,2,2,0,0）;Coefs = PP.coefs%3.两种差值函数图象比较figure（3）plot（X,Ln,-,color,1, 0, 0,LineWidth, 3）xlabel（深度x）;ylabel（温度T）;title（ Lagrange差值图）;holdonpause（3）fnplt（PP,b,3,0,28） %函数作图xlabel（深度x）;ylabel（温度T）;title（ 三次样条法差值图）;pause（3）plot（x,T,-,linewidth,3,markersize,10）运行结果:I_poly =（11156832321404503197312500000*t7942692507666801544174224450000*t6+29858157926737187739430003070000*t5 433725749060135044183186809635000*t4 + 2850751947133625442244870404829250*t38285950923799759686915051372346805*t2 + 8477820147617239232740036583612328*t+ 512657632901869980094703083834211328)/22484983899204823688364170343605760I_poly =Columns 1 through 6 0.0000-0.00000.0013-0.01930.1268-0.3685Columns 7 through 80.377022.8000Coefs =0.0022-0.0000-0.011522.8000-0.00920.01500.023022.8000-0.0114-0.0565-0.085022.80000.0297-0.2004-1.164020.6000-0.01530.2099-1.120013.90000.0017-0.0014-0.160511.7000-0.00170.0223-0.064011.1000图示如下：深度x深度x第六章 最小二乘拟合与最佳逼近 一、最小二乘拟合 加权最小二乘法逼近准则：min d (f, P ) = minPePPeP 艺 p f (xi=0 最小二乘逼近多项式P* (x)必须满足如下必要条件:n厂月、8 a088 a188an即满足法方程组:*t*a =*tfn a =那t*,1 *tf例题：下面是一处地质岩层断面高程(深度)的测量数据。水平距离km x ( ):00.201.002.103.505.006.807.509.0011.212.0亠心m,咼度():h1.641.581.681.841.580.860.390.310.390.770.86试利用最小二乘法求满足P*n2 0.5m=min|为(f (x )p (x )Ikn kzy(x) k=0即误差不超过0.5m的最低次数的拟合多项式P* (x)，写出该多项式的表达式;画出数据散点图n和该多项式曲线.程序：clc;clear all;close all;x=linspace(0,12,12);x=0 200 1000 2100 3500 5000 6800 7500 9000 11200 12000; h=1.64 1.58 1.68 1.84 1.58 0.86 0.39 0.31 0.39 0.77 0.86;plot(x,h,*,markersize,8) %画出给出数据的散点图xlabel(水平距离x);ylabel(高度h);title(岩层断面水平距离x和高度h的散点图);figure(1)%1.最小二乘法拟合n=3;P, S = polyfit(x, h, n)Pm = polyval(P, x)R = sqrt(sum(h - Pm).入2)% 误差t = linspace(x(1), x(end), 12);Poly = polyval(P, t);figure(2)plot(x, h, ro, t, Poly, LineWidth, 3, markersize, 8) set(gca,FontSize,18)legend(The Data, The Fitting Curve, 1)title(Curve Fitting by Least Square Approximation, fontsize, 18)0.20 2000 4000 6000 8000 10000 12000 水平距离XCurve Fitting by Least Square Approximation第七章 微积分的数值方法一、数值微分如果给定函数的关系式比较复杂或者未知，而仅仅知道在n +1个相异点X，k二0 1n k处的函数值，那么，我们可以利用函数的插值多项式的导数作为函数导数的近似值f (x)u L (x)=X l (x)f (x )nkkk=0因而有这里需要说明一点的是，尽管和的函数值可能相差不多，但是它们的导数有可能相差很大。二、数值积分数值积分方法最大的优点是将复杂的函数积分转化为相对简单的加、减、乘、除运算。对于给定的函数，如果 的函数关系式比较复杂或未知，而仅仅知道该函数在n +1个点处的函数值y，那么可以利用函数的插值多项式的积分作为函数在a, b上的定积分的近似 k值。1、Newton-Cotes 求积公式根据 Lagrange 插值多项式L (x)=X l (x)f (x )n k k k=0有I = fh f (x )dxfh L (x)dxl (x )dx f (x )kkA =f hl (x)dxak得 Newton-Cotes 求积公式fhf (x)dxfhL (x)dxaa n=艺f (x )f l (x)dxkkk =0a= LAf (x ) kk k=0特别地，当取插值节点为x, = a + kh, k = 0, 1,，nk,h - ah =n时有 两点公式（梯形公式）:fhf (x)dx fhL (x)dxa a n丄 Af (x ) kk k=0=亍:f (a )+ f (h )三点公式(Simpson公式)：anh a6fhf (x)dx fhL (x)dx+ f (h )a例题：下面是一处地质岩层断面上部边缘的深度测量数据。水平距离x00.201.002.103.505.006.807.509.0011.212.0(km ):深度y (m ):1.641.581.681.841.580.860.390.310.390.770.86表1试利用复化的梯形求积法求该组数据所在曲线与基准线(x轴)在范围o, 12内所围成图形面积.画出数据散点图和图形的示意图.试利用复化的Simpson求积法求该组数据所在曲线与基准线(x轴)在范围0, 9内所围成图形面积.画出数据散点图和图形的示意图.解答如下：程序：% 题目给出的数据XI = 0 0.20 1.00 2.10 3.50 5.00 6.80 7.50 9.00 11.2 12.00;YI = -1.64 1.58 1.68 1.84 1.58 0.86 0.39 0.31 0.39 0.77 0.86;n = length(XI) - 1;% 此时 n=10M = 60;x = linspace(XI(1), XI(end), M +1);%设置线性空间% 问题1.1T1 = 0; %初值%将图像分成若干个小区间，在每个区间内，求小梯形的面积，并累加起来，就是题目所要求的面积 for k = 2 : n + 1T1 = T1 + (YI(k) + YI(k-1)*(XI(k) - XI(k-1)/2;enddisp(复化的梯形求积法求得的图形面积：)T_Trape = T1%累加后的图形面积figure(1)%画出函数图像set(gca,FontSize, 20)patch(XI(1), XI, XI(end), 0, YI, 0,0.5 1 0.5)%颜色的填充holdonplot(XI, YI,o-, XI, 0*XI, k,LineWidth, 3,markersize, 8) title(复化的梯形求积法求围成图形面积)xlabel(水平距离x)ylabel(深度y)运行结果：复化的梯形求积法求得的图形面积：T_Trape =-11.6090图像：复化的梯形求积法求围成图形面积程序：函数 M 文件：function Ln = Lagrange_Fun_01(x, XI, YI) if isa(x, sym) = 1;n = length(XI) - 1;Ln = 0;Pn = sym(ones(n + 1, 1);for k = 1 : n + 1for i = 1 : n + 1 if i = kPn(k) = Pn(k)*(x - XI(i)/(XI(k) elsePn(k) = Pn(k);endendLn = Ln + Pn(k)*YI(k);endelsen = length(XI) - 1;L = ones(n + 1, length(x);Ln = zeros(size(x);for k = 1 : n + 1for i = 1 : n + 1if i = kL(k,:) = L(k,:).*(x - XI(i)./(XI(k) elseL(k,:) = L(k,:);endendLn = Ln + L(k, :).*YI(k);XI(i);XI(i);end end主程序：% 问题 1.2 和 1.3%问题3是问题2的延伸S2 = 0;%Simpson复化求积法，应用同样的面积累加，得到最终的图形面积for k = 1: 2: n - 1S2 = S2 + (XI(k+2) - XI(k)*(YI(k) + 4*YI(k+1) + YI(k+2)/6; X(k, :) = linspace(XI(k), XI(k+2), 20);Y(k, :) = Lagrange_Fun_01(X(k,:), XI(k : k+2), YI(k : k+2); enddisp( Simpson复化求积法求得的图形面积：)S_Simpson = S2 %Simpson复化求积法求得的面积figure(2) set(gca,FontSize, 20)holdonfor k = 1: 2 : n-1 patch(XI(k), X(k,:), XI(k +2), 0, Y(k, :), 0,0.5 1 0.5)plot(X(k, :), Y(k, :),LineWidth, 3)endplot(XI, YI,o, XI, 0*XI, k, LineWidth, 3, markersize, 8) title(复化的Simpson求积法求围成图形面积)xlabel(水平距离x)ylabel(深度y)运行结果：Simpson复化求积法求得的图形面积：S_Simpson =-11.9128图像：第八章微分方程（组）初值问题数值方法引言：微分方程（组）是科学研究和工程应用中最常用的数学模型之一。但是基绝大多数的常微分方程和常微分方程组得不到解析解和实际应用中往往只需要知道常微分方程（组）的解在某些点处的函数值这两点原因，绝大多数在实际中遇到的常微分方程和常微分方程组得不到“解析解”，我们可以采用下面将介绍的常微分方程（组）的初值问题的数值解法，就可以达到这一目的。Euler 方法最简单的数值解法1、Euler 法假设要求在点t = t + kh，k = 1,2,nk 0处初值问题（4）的解的近似值。首先对式（4）的两端积分，得对于该式的右边，如果被积函数用积分下限处的函数值代替被积函数作积分，则有y (t)= y (t )+J tk+1 f C, y (t)lzt1k ,tk-y (t )+ hft , y (t )kk k进而得到下式给出的递推算法一Euler方法y = y + hf (t , y ) k = 1,2,nk+1kk ky(t )= y002、Euler 方法的改进 改进的 Euler 方法 Euler 方法的精度高，其原因在于：1）在推导 Euler 方法时，我们是用待求解函数在一点处的变化率f （t , y （t ）kk 代替在区间上的平均变化率；而在推导改进的 Euler 方法时，我们是用待求解函数在两点处变化率的平均值(2)f (t, y (t )+ f (t , y (t )2k kk +1k +1代替在区间上的平均变化率； 显然，通常式（2）比式（1）更接近于在区间上的平均变化率。由此，人们受到启发：适当 地选取区间上函数多个点处的变化率，用它们加权平均值代替在区间上的平均变化率，近似 解的精度应更高。RungeKutta法：二阶的RungeKutta法三阶的RungeKutta法四阶的RungeKutta 法RungeKutta 法是按选取区间上函数变化率数目的多少和截断误差的阶数来区分的一系列 方法。例题：给定如下常微分方程 Cauchy 问题：dy = -ay 一 12 sin te-0.5/ + b, t e To, T* dtl fy(o)= y -0分别取参数，初值和步长：a = 1, b = 1 ，yo=-10.1 h 0.3。，在区间0 T 上分别利用Euler法，4阶R-K公式和mat lab的ode45求问题(1)的数值解；程序：clcclear alla=1;b=1;F = (t, y) -a*y-t2*exp(-0.5*t)*sin(t)+b;%欧拉法、matlab的ode45方法求解及画图y0 = -1;Tf = 2*pi;Tspan = 0, Tf;Options = odeset(RelTol, 1e-5, AbsTol, 1e-6);T, z = ode45(F, Tspan, y0,Options);h = 0.2;t1 = 0:h:Tf;Lt = length(t1);Y(1) = y0;for k = 2 : LtY(k) = Y(k-1) + h*F(t1(k-1), Y(k-1);endset(gca,Fontsize, 24)plot(T, z, t1, Y, t1, Y,o,LineWidth, 4, markersize, 5) %4阶R-K公式求解及画图tt = 0 : h : Tf;c = 1;YY(1) = y0;for k = 1: length(tt) -1;K1(k) = F(tt(k), YY(k);z(k) = YY(k) + 2h*K1(k);K2(k) = F(tt(k)+h/2, z(k);z(k) = YY(k) + 2h*K2(k);K3(k) = F(tt(k)+h/2, z(k);z(k) = YY(:,k) + h*K3(k);K4(k) = F(tt(k)+h, z(k);YY(k+1) = YY(k) + 6h*(K1(k) + 2*K2(k) + 2*K3(k) + K4(k);endhold onplot(tt, YY, LineWidth, 3)title(An Example of Applying The Fourth-Order Rong - Kutta Method) legend(Euler Curve, The Curve Obtained By Order 4 R-K,The Curve Obtained By ode45) xlabel(Time) ylabel(Orbit)%输出三种方法最后的数值解C=cell(2,3);C=Euler法,4 阶R-K法,ode45法；Y(end),YY(end),z(end)运行结果：最终得到的三种方法的数值解分别为：C =Euler法4阶 R-K 法ode4 5法 2.1863 2.0838 2.0029An Example of Applying The Fourth-Order Rong - Kutta Method参考文献：1、盖如栋，Leeture 1_引言、近似计算与误差分析，电子讲义，淮阴工学院，2014；2、盖如栋，数值分析实验，电子讲义，淮阴工学院，2013。