插值与数值积分PPT学习教案

1、会计学1 插值与数值积分插值与数值积分 计算机会“算”吗？靠得住吗 ？ 例：把例：把4开开n次方，再平方次方，再平方n次，结果是次，结果是4？存在误差？存在误差 ？ 英国著名数值分析学家英国著名数值分析学家 Higham (1998): Higham (1998): Can you count on computers?Can you count on computers? 精确计算：精确计算： 解析结果解析结果 (Analytical) 近似近似计计算：算： 数值数值结果结果 (Numerical ) ?4 2 2 n=55左右：结果变成左右：结果变成 1 计计算功效算功效= 计算工具计算工

2、具* 计算方法计算方法(算法算法) 浮点运算：舍入误差浮点运算：舍入误差 第1页/共42页 实验3的基本内容 3.3.数值积分的数值积分的梯形公式、辛普森公式和梯形公式、辛普森公式和高斯公式。高斯公式。 1.1.插值的基本原理；插值的基本原理； 三种插值方法：拉格朗日插三种插值方法：拉格朗日插 值，分段线值，分段线 性性 插值，三次样条插值。插值，三次样条插值。 2.2.插值的插值的 MATLAB 实现实现及插值的应用及插值的应用。 4.4.数值积分的数值积分的 MATLAB 实现实现及数值积分的应用及数值积分的应用。 第2页/共42页 什么是插值什么是插值(Interpolation)？从查

3、函数表说？从查函数表说 起起 查查 函函 数数 表表 xt dtex 2 2 2 1 )( x 012 1.0 0.8413 0.8438 0.8461 1.1 0.8643 0.8665 0.8686 1.2 0.8849 0.8869 0.8888 标准正态分布函数表标准正态分布函数表 求求 (1.114) (1.114)=0.8665 (0.8686 0.8665) 0.4=0.8673 插值插值 插值在图像处理插值在图像处理/数控加工数控加工/外观设计等领域有重要应用外观设计等领域有重要应用 第3页/共42页 插值的基本原理插值的基本原理 插值问题的提法插值问题的提法 已知已知 n+1

4、n+1个节点个节点, 1 , 0(),(njyx jj 其中其中 j x 互不相同，不妨设互不相同，不妨设), 10 bxxxa n 求任一插值点求任一插值点)( * j xx 处的插值处的插值. * y 0 x 1 x n x 0 y 1 y 节点可视为由节点可视为由 )(xgy 产生产生 , g表达式复杂表达式复杂, 甚至无表达甚至无表达 式式 * x * y 第4页/共42页 0 x 1 x n x 0 y 1 y 求解插值问题的基本思路求解插值问题的基本思路 构造一个构造一个( (相对简单的相对简单的) )函数函数),(xfy 通过全部节点通过全部节点, ,即即 ), 1 ,0()(n

5、jyxf jj 再用再用)(xf计算插值，计算插值， 即即 ).( * xfy * x * y 插值的插值的 基本原理基本原理 第5页/共42页 1.1.拉格朗日拉格朗日(Lagrange)(Lagrange)多项式插值多项式插值 1.0 1.0 插值多项式插值多项式 ) 1 ()( 01 1 1 axaxaxaxL n n n nn n n n n n n nn y y Y a a A xx xx X 0 0 1 1 00 , 1 1 在什么条件下）(0)det(X ), 1 , 0()(njyxL jjn )2(YXA 求 i a 三种插值三种插值 方法方法 有唯一解)2( 第6页/共42

6、页 1.1 1.1 拉格朗日插值多项拉格朗日插值多项 式式 ni xxxxxxxx xxxxxxxx xl niiiiii nii i 1 , 0, )()()( )()()( )( 110 110 )3()()( 0 xlyxL i n i in jjnji yxL ji ji xl )( ,0 , 1 )(则若 又又（2）有唯一解，故有唯一解，故（3）与与（1）相同相同 。 基函数基函数 ( ) i lx ) 1 ()( 01 1 1 axaxaxaxL n n n nn ) 2(YXA 三种插值三种插值 方法方法 第7页/共42页 ),(),( )!1( )( )()()( 0 )1(

7、baxx n g xLxgxR n j j n nn 1 )1( )( n n Mg 减小（粗略地看）如何使误差)(xRn 平缓g j xx 接近 n j j n n xx n M xR 0 1 )!1( )( 三种插值三种插值 方法方法 1.2 1.2 误差估计误差估计 增加n 第8页/共42页 1.3 1.3 拉格朗日插值多项式的振荡拉格朗日插值多项式的振荡 ?)(?)(xRxLn nn 55, 1 1 )( 2 x x xg 63. 363. 3),()(lim xxgxLn n Runge现现 象象 取n=2,4,6,8,10, 计算Ln(x), 画出 图形 -505 -1.5 -1

8、-0.5 0 0.5 1 1.5 2 y=1/(1+x2) n=2 n=4 n=6 n=8 n=10 三种插值三种插值 方法方法 Runge.m 第9页/共42页 2.2.分段线性插值分段线性插值 xjxj-1xj+1x0 xn 其它,0 , , )( )()( 1 1 1 1 1 1 0 jj jj j jj jj j j n j jjn xxx xx xx xxx xx xx xl xlyxI 计算量与计算量与n n无关无关; ; n n越大，误差越小越大，误差越小 . . nn n xxxxgxI 0 ),()(lim 三种插值方法三种插值方法 第10页/共42页 机翼下轮廓线 3. 3

9、. 三次样条插值三次样条插值 样条函数的由来样条函数的由来 飞机、船体、汽车外形等的放样（设计飞机、船体、汽车外形等的放样（设计 ） 细木条：样条细木条：样条 第11页/共42页 3. 3. 三次样条插值三次样条插值 , 1,),()( 1 nixxxxsxS iii ,)()3 ), 1 ,0()()2 ), 1()()1 0 2 23 n ii iiiii xxCxS niyxS nidxcxbxaxs 数学样条（数学样条（spline） iiii dcba n , 4 个待定系数 3 ） ) 1, 1()()( )()(),()( 1 11 nixsxs xsxsxsxs iiii ii

10、iiiiii 3） 2），3）共 4n-2个方程 三种插值方法三种插值方法 第12页/共42页 自然边界条件）(0)()()4 0 n xSxS )(,)4)3)2xSdcba iiii 三次样条插值确定三次样条插值确定4 4n n个系数需增加个系数需增加 2 2个条件个条件 思考1）自然边界条件的几何意义是什么 ？ 2）样条插值为什么普遍用3次多项式 ，而不是2或4次？ 三次样条插值三次样条插值 ).()(limxgxS n 第13页/共42页 三种插值方法小结三种插值方法小结 拉格朗日插值（高次多项式插值）：拉格朗日插值（高次多项式插值）： 曲线光滑；误差估计有表达式；收敛性不能保证。曲线

11、光滑；误差估计有表达式；收敛性不能保证。 用于理论分析，实际意义不大用于理论分析，实际意义不大。 分段线性和三次样条插值（低次多项式插值）：分段线性和三次样条插值（低次多项式插值）： 曲线不光滑（三次样条插值已大有改进）；误差估曲线不光滑（三次样条插值已大有改进）；误差估 计较难（对三次样条插值）；收敛性有保证。计较难（对三次样条插值）；收敛性有保证。 简单实用，应用广泛简单实用，应用广泛。 其他：其他：Hermite插值、分段三次插值、二维插值等插值、分段三次插值、二维插值等 根据需要，根据需要，各取所需各取所需。 第14页/共42页 1. 1. 拉格朗日插值拉格朗日插值: :自编程序自编程

12、序, ,如名为如名为 lagr.m 的的M文件，文件， 第一行为第一行为 function y=lagr(x0,y0,x) 输入输入: :节点节点x0,y0, 插值点插值点x ( (均为均为数组，长度自定义数组，长度自定义) )）；）； 输出输出: :插值插值y ( (与与x同长度数组同长度数组) )）。）。 应用时输入应用时输入x0,y0,x后后, ,运行运行 y=lagr(x0,y0,x) 2. 2. 分段线性插值分段线性插值: :已有程序已有程序 y=interp1(x0,y0,x) y=interp1(x0,y0,x,linear) 3. 3. 三次样条插值三次样条插值: :已有程序已

13、有程序 y=interp1(x0,y0,x,spline) 或或 y=spline(x0,y0,x) 用MATLAB作插值计算 注：注：MATLAB有样条工具箱（有样条工具箱（Spline Toolbox） 第15页/共42页 用MATLAB作插值计算 55, 1 1 )( 2 x x xg 为例，作三种插值的比较为例，作三种插值的比较以以 0 1.0000 1.0000 1.0000 1.0000 0.5000 0.8000 0.8434 0.7500 0.8205 1.0000 0.5000 0.5000 0.5000 0.5000 1.5000 0.3077 0.2353 0.3500

14、0.2973 2.0000 0.2000 0.2000 0.2000 0.2000 2.5000 0.1379 0.2538 0.1500 0.1401 3.0000 0.1000 0.1000 0.1000 0.1000 3.5000 0.0755 -0.2262 0.0794 0.0745 4.0000 0.0588 0.0588 0.0588 0.0588 4.5000 0.0471 1.5787 0.0486 0.0484 5.0000 0.0385 0.0385 0.0385 0.0385 x y y1 y2 y3 用用n=11个节点，个节点，m=21个插值点，三种方法作插值，画图。

15、个插值点，三种方法作插值，画图。 chazhi1 第16页/共42页 插值的应用 加工时需要加工时需要x每改变每改变0.05时的时的y值值 chazhi2 图1 零件的轮廓线 （x间隔0.2） 表1 x间隔0.2的加工坐标x,y（图1右半部的数据） 数控机床加工零件数控机床加工零件 0.0,5.00 0.2,4.710.4,4.31 0.6,3.68 0.8,3.05 1.0,2.50 1.2,2.051.4,1.69 1.6,1.40 1.8,1.18 2.0,1.00 2.2,0.862.4,0.74 2.6,0.64 模型模型 将图1逆时针方向转90度，轮廓线上下对称，只需对上半部计算一

16、个函数在插值点的值。 图2 逆时针方向转90度的结果 -5-4-3-2-1012345 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 u v 令v=x, u= -y 第17页/共42页 为什么要作数值积分为什么要作数值积分 许多函数许多函数“积不出来积不出来”, ,只能用数值方法，如只能用数值方法，如 dx x x dxe b a b a x sin , 2 2 积分是重要的数学工具，是微分方程、概率积分是重要的数学工具，是微分方程、概率 论等的基础；在实际问题中有直接应用。论等的基础；在实际问题中有直接应用。 对于用离散数据或者图形表示的函数对于用离散数据或者图形表示的函数

17、，(可以先做插值然后积分；或者直接利用点做数值积分) 计算积分只有求助于数值方法。计算积分只有求助于数值方法。 数值积分数值积分 第18页/共42页 n ab fIIdxxfI n k knn n b a )(,lim)( 1 数数 值值 积积 分分 的的 基基 本本 思思 路路 回回 忆忆 定定 积积 分分 的的 定定 义义 各种数值积分方法研究的是各种数值积分方法研究的是 k ),(ba如何取值，区间如何取值，区间如何划分，如何划分， 使得既能保证一定精度，计算量又小。使得既能保证一定精度，计算量又小。 n n充分大时充分大时I In n就是就是I I的数值积分的数值积分 （计算功效：算得

18、准，算得快）（计算功效：算得准，算得快） 第19页/共42页 1.1.从矩形公式到梯形公式数值积分数值积分 y y=f(x) xb a o )1( 1 0 n k kn fhL )(, , 10 kk nk xff n ab h bxxxxa )2( 1 n k kn fhR nn RL,平均，得到 梯形公式) 3()( 2 0 1 1 n n k kn ff h fhT xk+ 1 x k xk-1 fk 第20页/共42页 2.2.辛普森辛普森(Simpson)(Simpson)公式公式 （抛物线公式）（抛物线公式） 梯形公式相当于用分段线性插值函数分段线性插值函数代替)(xf 每段要用相

19、邻两小区间两小区间 端点的三个函数值端点的三个函数值 抛物线抛物线 公式公式 提高精度提高精度分段二次插值函数分段二次插值函数 2221212222 (,),(,),(,) 0,1, ,1 kkkkkk xfxfxf km 数值积分数值积分 y y=f(x) xb a o x 2k f2k x2k+1x2k+2 f2k+ 1 f2k+2 区间数必须为偶数区间数必须为偶数mn2 第21页/共42页 ) 4( 2 ),24( 3 1 1 2 1 0 1220 m ab hffff h S m k k m k kmm 对对k求和求和(共共m段段) )，得（复合），得（复合）辛普森公式辛普森公式 ：

20、)4( 3 )( 22122 22 2 kkk x x k fff h dxxs k k 二次插值函数sk(x) 构造用),(),(),( 2222121222kkkkkk fxfxfx 2.2.辛普森辛普森(Simpson)(Simpson)公式（抛物线公式）公式（抛物线公式） 第22页/共42页 b a nnn TdxxfTITfR)(),( 梯形公式在每小段上是用梯形公式在每小段上是用线性插值函数线性插值函数T T( (x) )代替代替 f( (x) ) ),(,),)( 2 )( )()( 11 kkkkk k xxxxxxx f xTxf 梯形公式梯形公式 的误差估计的误差估计 )(

21、 2 0 1 1 n n k kn ff h fhT b a dxxf)( )( 12 )( 2 )( )()( 3 1 11 k x x kk k x x f h dxxxxx f dxxTxf k k k k 因为：因为：(x-xk)(x-xk+1)在在(xk,xk+1)不变号，所以：不变号，所以： 第23页/共42页 )5()( 12 |),(| 2 2 abM h TfR n 梯形公式梯形公式Tn的的 误差是误差是h2阶的阶的 ),(, )(max 2 baxxfM 估计估计 h ab n 因为 1 0 3 )( 12 |),(| n k kn f h TfR 梯形公梯形公 式的误差式

22、的误差 )( )( ( 12 1 )( 12 1 2 afbfdxxf h TI b a n 1 0 3 )( 12 )( n k k b a nn f h TdxxfTI ) 5()()( 12 2 afbf h TI n 第24页/共42页 同理可得：同理可得： ) 6()( 180 | ),(| 4 4 abM h SfR n 其中其中 ),(,)(max )4( 4 baxxfM 辛普森公式辛普森公式Sn的误差是的误差是h4阶的阶的。 辛普森公式的误差估辛普森公式的误差估 计计 第25页/共42页 梯形公式和辛普森公式的收敛性 若对若对I某个数值积分某个数值积分In有有c h II p

23、 n n lim（非零常数）（非零常数） 则称则称 In是是 p 阶收敛的阶收敛的。 梯形公式梯形公式 2 2 阶收敛，辛普森公式阶收敛，辛普森公式 4 4 阶收敛。阶收敛。 c=0: 至少至少p阶收敛（超阶收敛（超p阶收敛）阶收敛） 第26页/共42页 积分步长的自动选取积分步长的自动选取 选定数值积分公式后，如何确定步长选定数值积分公式后，如何确定步长h以满足给定的误差以满足给定的误差 )()( 12 2 afbf h TI n 梯形公式 )( 4 1 2nn TITI nn TT 2 用二分法只要用二分法只要 其中其中fk+1/2是原是原分点分点xk,xk+1的中点的中点（记记xk+1/

24、2）的函数值的函数值 1 0 2 12 22 n k k n n f hT T 且且T T2n 2n可在 可在T Tn n 基础上计算基础上计算 )( 3 1 22nnn TTTI n TI 2 )2( 2 nn h h 第27页/共42页 高斯高斯(Gauss)(Gauss)求积公式求积公式 矩形公式矩形公式(1)、(2) 梯形公式梯形公式(3) 辛普森公式辛普森公式(4) A Ak k是与是与f f无关的常数无关的常数 代数代数 精度精度 设设 ,)( k xxf 用用(7)计算计算 ,)( b a dxxfI 若对于若对于mk,1 ,0 都有都有, IIn 而当而当, 1IImk n 则

25、则称称In的代数精度为的代数精度为m. )7()( 1 n k kkn xfAI Newton-Cotes方法方法 第28页/共42页 梯形公式的代数精度（考察梯形公式的代数精度（考察T1） k=1 f(x)=x 2 22 ab xdxI b a 2 )( )()( 2 1 baab bfaf h T 3 33 2 ab dxxI b a 2 )( 22 1 baab T k=2 f(x)=x2 IT 1 IT 1 梯形公式的代数精度为梯形公式的代数精度为1辛普森公式的代数精度为辛普森公式的代数精度为3 第29页/共42页 高斯公式的思路高斯公式的思路 取消对节点的限制，按照代数精度最大取消对

26、节点的限制，按照代数精度最大 的原则，同时确定节点的原则，同时确定节点xk和系数和系数Ak 构造求积公式构造求积公式 )()( 22112 xfAxfAG 对于对于 1 1 )(dxxfI 使使G G2 2的代数精度为的代数精度为3 3 32, , 1)(xxxxf )()()( 2211 1 1 xfAxfAdxxf 确定确定 2121 ,AAxx 第30页/共42页 0 3/2 0 2 3 2 2 3 1 1 2 22 2 11 2211 21 xAxA xAxA xAxA AA 将将f(x)f(x)代入计算得代入计算得 1,3/1,3/1 2121 AAxx )3/1 ()3/1( 2

27、ffG 用用n个节点，个节点，Gn的代数精度可达的代数精度可达2n-1, 但是需解但是需解 复杂的非线性方程组，实用价值不大。复杂的非线性方程组，实用价值不大。 第31页/共42页 常 用 的 高 斯 公 式 将将( (a,ba,b) )分小，把小区间变换为分小，把小区间变换为(-1(-1，1), 1), 再用再用G G2 2 m k kk b a zfzf h dxxf 1 )2()1 ( )()( 2 )( 32 2 , 32 2 1) 2(1) 1 ( hxx z hxx z kk k kk k mkkhaxmabh k , 1, 0,/ )( 代数精度为代数精度为3节点加密时，原计算信

28、息无法利用节点加密时，原计算信息无法利用 第32页/共42页 思路思路：将积分区间分小，在小区间上用：将积分区间分小，在小区间上用n不太不太 大大 的的 。而在节点加密一倍时能够利用原节点的。而在节点加密一倍时能够利用原节点的 函数值，可以把区间的端点作为固定节点。函数值，可以把区间的端点作为固定节点。 改进的高斯公式 n G )()()( 1 2 1 bfAxfAafAG n n k kkn Gauss-Lobatto求积公式求积公式 其中a, b为小区间的端点， nn AAxx, 112 为2n-2个参数， 代数精度可达到代数精度可达到2n-3 注意：实际计算中一般采用自适应方法确定步长注

29、意：实际计算中一般采用自适应方法确定步长 第33页/共42页 用用MATLAB 作数值积分作数值积分 1 0 n k kn fhL n k kn fhR 1 矩形矩形 公式公式 Sum(x) 输入数组x(即fk),输出x的和(数) cumsum(x)输入数组x,输出x的依次累加和(数组) 梯形梯形 公式公式 )( 2 0 1 1 n n k kn ff h fhT trapz(x) 输入数组x,输出按梯形公式x的积分(单位步长) trapz(x,y) 输入同长度数组 x,y,输出按梯形公式 y对x的积分(步长不一定相等) 第34页/共42页 用用MATLAB 作数值积分作数值积分 m ab h

30、ffff h S m k k m k kmn 2 ),24( 3 1 1 2 1 0 1220 辛普森公式辛普森公式 quad(fun,a,b,tol,trace) I,fn=quad() 用自适应辛普森公式计算 tol为绝对误差，缺省时为10-6 Gauss-Lobatto公公 式式 )()()( 1 2 1 bfAxfAafAG n n k kkn quadl(fun,a,b,tol,tr ace) I,fn=quadl() 用自适应Gauss-Lobatto公式计 算 tol为绝对误差，缺省时为10- 6 注意：注意：fun.m中应以自变量为矩阵的形式输入中应以自变量为矩阵的形式输入(点

31、运算点运算 ) 第35页/共42页 矩形域上计算二重积分的命令：矩形域上计算二重积分的命令： dblquad(fun,xmin,xmax,ymin,ymax,tol) 广义积分、二重和三重积分 长方体上计算三重积分的命令：长方体上计算三重积分的命令： triplequad(fun,xmin,xmax,ymin,ymax, zmin,zmax,tol) 注：注：fun是被积函数，本身可以有自己的参数是被积函数，本身可以有自己的参数 广义积分：广义积分： 通过分析和控制误差，转换成普通积分通过分析和控制误差，转换成普通积分 quadv(fun,a,b,tol,trace) 向量值积分：向量值积分： 第36页/共42页 用用MATLAB 作数值积分作数值积分 例例. 计算计算 4 0 1 1s i n d x x 1 1）矩形公式和梯