1 大学数学实验 Experiments in Mathematics 实验3 插值与数值积分 清华大学数学科学系 实验3的基本内容 3.数值积分的梯形公式、辛普森公式和高斯公式。 1.插值的基本原理； 三种插值方法：拉格朗日插值，分段线性 插值，三次样条插值。 2.插值的 MATLAB 实现及插值的应用。 4.数值积分的 MATLAB 实现及数值积分的应用。 什么是插值？ 从查函数表说起 查函数表 ∫ ∞? ? =Φ x t 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 插 值 插值的基本原理 插值问题的提法 已知 n+1个节点 ,,1,0(),( njyx jj L= 其中 j x 互不相同，不妨设 ), 10 bxxxa n =<<<= L 求任一插值点 )( * j xx ≠ 处的插值 . * y ? ? ? ? ? 0 x 1 x n x 0 y 1 y 节点可视为由 )( xgy = 产生 , g表达式复杂 , 甚至无表达式 ? * x * y ? ? ? ? ? 0 x 1 x n x 0 y 1 y 求解插值问题的基本思路 构造一个(相对简单的)函数 ),(xfy = 通过全部节点,即 ),1,0()( njyxf jj L== 再用 )(xf 计算插值，即 ).( ** xfy = ? * x * y 插值的 基本原理 1.拉格朗日(Lagrange)多项式插值 1.0 插值多项式 )1()( 01 1 1 axaxaxaxL n n n nn ++++= ? ? L ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? = ? ? n n n n n n nn y y Y a a A xx xx X MM L L L 0 0 1 1 00 ,, 1 1 在什么条件下）(0)det( ≠XQ ),1,0()( njyxL jjn L== )2(YXA = 求 i a 三种插值 方法 有唯一解)2(∴ 2 1.1 拉格朗日插值多项式 ni xxxxxxxx xxxxxxxx xl niiiiii nii i L LL LL 1,0, )())(()( )())(()( )( 110 110 = ???? ???? = +? +? )3()()( 0 xlyxL i n i in ∑ = = jjnji yxL ji ji xl =∴ ? ? ? ≠ = = )( ,0 ,1 )(Q 又（ 2）有唯一解，故（ 3）与（ 1）相同。 基函数 () i lx )1()( 01 1 1 axaxaxaxL n n n nn ++++= ? ? L )2(YXA= 三种插值 方法 ),(),( )!1( )( )()()( 0 )1( baxx n g xLxgxR n j j n nn ∈? + =?= ∏ = + ξ ξ 1 )1( )( + + ≤ n n Mg ξ 减小（粗略地看）如何使误差 )(xR n 平缓g j xx 接近 ∏ = + ? + ≤ n j j n n xx n M xR 0 1 )!1( )( 三种插值 方法 1.2 误差估计 增加n 1.3 拉格朗日插值多项式的振荡 ?)(?)( ↓??↑ xRxLn nn 55, 1 1 )( 2 ≤≤? + = x x xg 63.363.3),()(lim ≤≤?= ∞→ xxgxL n n Runge现象 取 n=2,4,6,8,10,计 算 L n (x), 画出图形 -5 0 5 -1.5 -1 -0.5 0 0.5 1 1.5 2 y=1/(1+x 2 ) n=2 n=4 n=6 n=8 n=10 三种插值 方法 Matlab.lnk Runge.m 2.分段线性插值 ? ? ? ? ? ? x j x j-1 x j+1 x 0 x n ? ? ? ? ? ? ? ? ? ? ? ≤≤ ? ? ≤≤ ? ? = = + + + ? ? ? = ∑ 其它,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越大，误差越小. nn n xxxxgxI ≤≤= ∞→ 0 ),()(lim 三种插值 方法 机翼下轮廓 线 3. 三次样条插值 样条函数的由来 飞机、船体、汽车外形等的放样（设计） 细木条：样条 3. 三次样条插值 },1],,[),({)( 1 nixxxxsxS iii L=∈= ? ],[)()3 ),1,0()()2 ),1()()1 0 2 23 n ii iiiii xxCxS niyxS nidxcxbxaxs ∈ == =+++= L L 数学样条（ spline） iiii dcba n ,,, 4 个待定系数 3） )1,1()()( )()(),()( 1 11 ?=′′=′′ ′=′= + ++ nixsxs xsxsxsxs iiii iiiiiiii 3 ’ ） 2）， 3 ’ ）共 4n-2个方程 三种插值 方法 3 自然边界条件）(0)()()4 0 =′′=′′ n xSxS )(,,,)4)3)2 xSdcba iiii ?? 三次样条插值确定4n个系数需增加 2个条件 思考 1）自然边界条件的几何意义是什么？ 2）样条插值为什么普遍用3次多项式， 而不是2或4次？ 三次样条 插值 ).()(lim xgxS n = ∞→ 三种插值方法小结 ? 拉格朗日插值（高次多项式插值）： 曲线光滑；误差估计有表达式；收敛性 不能保证（振荡现象）。 用于理论分析，实际意义不大 。 ? 分段线性和三次样条插值（低次多项式插值）： 曲线不光滑（三次样条插值已大有改进）；误差估 计较难（对三次样条插值）；收敛性有保证。 简单实用，应用广泛 。 1. 拉格朗日插值:自编程序,如名为 lagr.m 的 M文件， 第一行为 function y=lagr(x0,y0,x) 输入:节点 x0,y0, 插值点 x (均为数组，长度自定义)）； 输出:插值 y (与 x同长度数组)）。 应用时输入 x0,y0,x后,运行 y=lagr(x0,y0,x) 2. 分段线性插值:已有程序 y=interp1(x0,y0,x) 3. 三次样条插值:已有程序 y=interp1(x0,y0,x,’spline’) 或 y=spline(x0,y0,x) 用 MATLAB作插值计算 注： lagr.m 程序 可参考课本； MATLAB有样条工具箱（ Spline Toolbox） 用 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 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个 插值点，三 种方法作插 值，画图。 Matlab.lnk chazhi1 插值的应用 加工时需要 x每 改变 0.05时的 y值 MATLAB 5.3.lnk chazhi2 图1 零件的轮廓线 （ x间隔 0.2） 表 1 x间隔 0.2的加工坐标 x,y（图 1右半部的数据） 数控机床加工零件 ………2.6,0.642.4,0.742.2,0.862.0,1.00 1.8,1.181.6,1.401.4,1.691.2,2.051.0,2.50 0.8,3.050.6,3.680.4,4.310.2,4.710.0,5.00 模型 将图 1逆时针方向转 90度，轮廓线上下 对称，只需对上半部计算一个函数在插值点 的值。 图 2 逆时针方向转 90度的结果 -5 -4 -3 -2 -1 0 1 2 3 4 5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 u v 令 v=x, u=-y 为什么要作数值积分 ? 许多函数“积不出来”,只能用数值方法，如 dx x x dxe b a b a x ∫∫ ? sin , 2 2 ? 积分是重要的数学工具，是微分方程、概率 论等的基础；在实际问题中有直接应用。 ? 对于用离散数据或者图形表示的函数 ， 计算积分只有求助于数值方法。 数值 积分 4 n ab fIIdxxfI n k knn n b a ? === ∑ ∫ = ∞→ )(,lim)( 1 ξ 数值积分的基本思路 回忆定积分的定义 各种数值积分方法研究的是 k ξ ),( ba如何取值，区间 如何划分， 使得既能保证一定精度，计算量又小。 n充分大时 I n 就是 I的数值积分 1.从矩形公式到梯形公式数值积分 y y=f(x) xb a o )1( 1 0 ∑ ? = = n k kn fhL )(, , 10 kk nk xff n ab h bxxxxa = ? = =<<<= LL )2( 1 ∑ = = n k kn fhR nn RL , 平均，得到 梯形公式 )3()( 2 0 1 1 n n k kn ff h fhT ++= ∑ ? = x k+1 x k x k-1 f k 2.辛普森(Simpson)公式 （抛物线公式） 梯形公式相当于用 分段线性插值函数 代替 )(xf 每段要用相邻 两小区间 端点的三个函数值 抛物线 公式 提高精度 分段二次插值函数 2 2 21 21 22 22 (,),( , ),( , ) 0, 1, , 1 kk k k k k xf x f x f km ++ ++ =?L 数值积分 y y=f(x) xb a o x 2k f 2k x 2k+1 x 2k+2 f 2k+1 f 2k+2 区间数必须为偶数 mn 2= )4( 2 ),24( 3 1 1 2 1 0 1220 m ab hffff h S m k k m k kmm ? =+++= ∑∑ ? = ? = + 对 k求和 (共 m段)，得 辛普森公式 ： )4( 3 )( 22122 22 2 ++ ++= ∫ + kkk x x k fff h dxxs k k 二次插值函数 s k (x)构造用 ),(),,(),,( 2222121222 ++++ kkkkkk fxfxfx 2.辛普森(Simpson)公式（抛物线公式） ∫ ?= b a nn TdxxfTfR )(),( 梯形公式在每小段上是用 线性插值函数 T(x)代替 f(x) ),(,),)(( 2 )( )()( 11 ++ ∈?? ′′ += kkkkk k xxxxxxx f xTxf η η 梯形公式 的误差估计 )( 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-x k )(x-x k+1 )在 (x k ,x k+1 )不变号，所以： )5()( 12 |),(| 2 2 abM h TfR n ?≤ 即 梯形公式 T n 的 误差是 h 2 阶的 ),(,)(max 2 baxxfM ∈′′= 估计 h ab n ? =因为 ∑ ? = ′′≤ 1 0 3 )( 12 |),(| n k kn f h TfR ξ 梯形公式 的误差 ))(')('( 12 1 )( 12 1 )( 12 )( '' 2 1 0 3 afbfdxxf h TI f h TdxxfTI b a n n k k b a nn ??=?→ ? ′′?=?=? ∫ ∑ ∫ ? = ξ 5 同理可得： )6()( 180 |),(| 4 4 abM h SfR n ?≤ 其中 ),(,)(max )4( 4 baxxfM ∈= 即 辛普森公式 S n 的误差是 h 4 阶的 。 辛普森公式的误差估计 梯形公式和辛普森公式的收敛性 若对 I某个数值积分 I n 有 c h II p n n = ? ∞→ lim （非零常数） 则称 I n 是 p 阶收敛的 。 梯形公式 2 阶收敛，辛普森公式 4 阶收敛。 积分步长的自动选取 选定数值积分公式后，如何确定步长 h以满足给定的误差 ε ))()(( 12 '' 2 afbf h TI n ??≈?梯形公式 )( 4 1 2 nn TITI ?≈? ε≤? nn TT 2 用二分法只要 其中 f k+1/2 是原分点 x k ,x k+1 的中点 （ 记 x k+1/2 ） 的函数值 ∑ ? = + += 1 0 2 12 22 n k k n n f hT T 且 T 2n 可在 T n 基础上计算 )( 3 1 22 nnn TTTI ?≈? ε≤? n TI 2 )2( 2 nn h h →→ 高斯(Gauss)求积公式 矩形公式 (1)、 (2) 梯形公式 (3) 辛普森公式 (4) A k 是与 f无关的常数 代数 精度 设 ,)( k xxf = 用 (7)计算 ,)( ∫ = b a dxxfI 若对于 mk L,1,0= 都有 ,II n = 而当 ,,1 IImk n ≠+= 则 称 I n 的代数精度为 m. )7()( 1 ∑ = = n k kkn xfAI Newton -Cotes 方法 梯形公式的代数精度（考察 T 1 ） 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)=x 2 IT = 1 IT ≠ 1 梯形公式的代数精度为 1 辛普森公式的代数精度为 3 高斯公式的思路 取消对节点的限制，按照代数精度最大 的原则，同时确定节点 x k 和系数 A k 构造求积公式 )()( 22112 xfAxfAG += 对于 ∫ ? = 1 1 )( dxxfI 使G 2 的代数精度为3 32 ,,,1)( xxxxf = )()()( 2211 1 1 xfAxfAdxxf += ∫ ? 确定 2121 ,,, AAxx 6 0 3/2 0 2 3 2 2 3 1 1 2 22 2 11 2211 21 =+ =+ =+ =+ xAxA xAxA xAxA AA 将f(x)代入计算得 1,3/1,3/1 2121 ===?= AAxx )3/1()3/1( 2 ffG +?= 用 n个节点， G n 的代数精度可达 2n-1, 但是需解 复杂的非线性方程组，实用价值不大。 常用的高斯公式 将( a,b)分小，把小区间变换为(-1，1), 再用 G 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 L,1,0,,/)( =+=?= 代数精度为 3 思路：将积分区间分小，在小区间上用 n不太 大 的 。而在节点加密一倍时能够利用原节点的函 数值，需要把区间的端点作为固定节点。可以改 进高斯公式 改进的高斯公式 n G )()()( 1 2 1 bfAxfAafAG n n k kkn ++= ∑ ? = Gauss-Lobatto求积公式 其中 a, b为小区间的端点， nn AAxx LL ,,,, 112 ? 为 2n-2个参数， 代数精度可达到 2n-3 用 MATLAB 作数值积分 ∑ ? = = 1 0 n k kn fhL ∑ = = n k kn fhR 1 矩形 公式 Sum(x) 输入数组 x(即 f k ),输出 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的积分(步长不一定相等) 用 MATLAB 作数值积分 m ab hffff h S m k k m k kmn 2 ),24( 3 1 1 2 1 0 1220 ? =+++= ∑∑ ? = ? = + 辛普森公式 quad(@fun,a,b,tol) 用自适应辛普森公式计算 tol为绝对误差，缺省时为 10 -6 Gauss-Lobatto公式 )()()( 1 2 1 bfAxfAafAG n n k kkn ++= ∑ ? = quadl(@fun,a,b,tol) 用自适应 Gauss-Lobatto公式计算 tol为绝对误差，缺省时为10 -6 矩形域上计算二重积分的命令： dblquad(@fun,xmin,xmax,ymin,ymax,tol) fun是被积函数 二重和三重积分 长方体上计算三重积分的命令： triplequad(@fun,xmin,xmax,ymin,ymax, zmin,zmax,tol) fun是被积函数 7 用 MATLAB 作数值积分 例 . 计算 4 0 1 1sin dx x π ? ∫ 1）矩形公式和梯形公式 将(0, π /4)100等分 2）辛普森公式和 Gauss-Lobatto公式 精确、方便 无法计算用数值给出的函数的积分 MATLAB 5.3.lnk Jifen1.m 数值积分的应用 实例 人造卫星轨道长度 )20( sin,cos π≤≤ == t tbytax 决定由 短半轴 长半轴 rss b a ,, ~ ~ 21 dttbtadtyxL ∫∫ +=+= 2 0 2222 2 0 22 cossin44 ππ && 轨道长度 y x o 近地点 s 1 =439km,远地点 s 2 = 2384km s 1s 2 地球半径 r=6371km r 需要作数值积分 s 1 =439km, s 2 = 2384km, r=6371km y x o s 1s 2 r s 1 s 2 y x o r a c b 决定由短半轴长半轴 rssba ,,,~,~ 21 数值积分实例 人造卫星轨道长度 dttbtaL ∫ += 2 0 2222 cossin4 π 12 22 ssra ++= 7782.5 2 12 = + += ss ra 1 srac ??=焦距 2 12 ss c ? = 7721.5 22 =?= cab dttbtaL ∫ += 2 0 2222 cossin4 π 用梯形公式和辛普森公式计算 只将区间 5等分，梯形公式就给出很好的结果 轨道长度 L=4.8707×10 4 千米 数值积分实例 人造卫星轨道长度 MATLAB 5.3.lnk Jifen2.m 布置实验 目的 1、掌握用 MATLAB计算拉格朗日、分段线性、三次样 条三 种插值的方法，改变节点的数目，对三种插值 结果进行初步分析。 2、掌握用 MATLAB及梯形公式、辛普森公式计算数值 积分。 3、 通过实例学习用插值和数值积分解决实际问题。 内容 课上布置，或参见网络学堂