1、大学数学实验大学数学实验Mathematical Experiments 实验实验3 3 插值与数值积分插值与数值积分计算机会“算”吗？靠得住吗？例：把例：把4开开n次方，再平方次方，再平方n次，结果是次，结果是4？存在误差？存在误差？英国著名数值分析学家英国著名数值分析学家 HighamHigham(1998):(1998):Can you count on computers?Can you count on computers?精确计算：精确计算：解析结果解析结果(Analytical)近似计算：近似计算：数值结果数值结果(Numerical)?422n=55左右：结果变成左右：结果变成
2、1计算功效计算功效=计算工具计算工具*计算方法计算方法(算法算法)浮点运算：舍入误差浮点运算：舍入误差实验3的基本内容3.3.数值积分的数值积分的梯形公式、辛普森公式和梯形公式、辛普森公式和高斯公式。高斯公式。1.1.插值的基本原理；插值的基本原理；三种插值方法：拉格朗日插三种插值方法：拉格朗日插 值，分段线性值，分段线性 插值，三次样条插值。插值，三次样条插值。2.2.插值的插值的 MATLAB 实现实现及插值的应用及插值的应用。4.4.数值积分的数值积分的 MATLAB 实现实现及数值积分的应用及数值积分的应用。什么是插值什么是插值(Interpolation)？从查函数表说起？从查函数表
3、说起查查 函函 数数 表表xtdtex2221)(x0121.0 0.8413 0.8438 0.84611.1 0.8643 0.8665 0.86861.2 0.8849 0.8869 0.8888标准正态分布函数表标准正态分布函数表求求 (1.114)(1.114)=0.8665 (0.8686 0.8665)0.4=0.8673插值插值插值在图像处理插值在图像处理/数控加工数控加工/外观设计等领域有重要应用外观设计等领域有重要应用插值的基本原理插值的基本原理插值问题的提法插值问题的提法已知已知 n+1n+1个节点个节点,1,0(),(njyxjj其中其中jx互不相同，不妨设互不相同，不
4、妨设),10bxxxan求任一插值点求任一插值点)(*jxx 处的插值处的插值.*y0 x1xnx0y1y节点可视为由节点可视为由)(xgy 产生产生,g表达式复杂表达式复杂,甚至无表达式甚至无表达式*x*y0 x1xnx0y1y求解插值问题的基本思路求解插值问题的基本思路构造一个构造一个(相对简单的相对简单的)函数函数),(xfy 通过全部节点通过全部节点,即即),1,0()(njyxfjj再用再用)(xf计算插值，即计算插值，即).(*xfy*x*y插值的插值的基本原理基本原理1.1.拉格朗日拉格朗日(Lagrange)(Lagrange)多项式插值多项式插值1.0 1.0 插值多项式插值
5、多项式)1()(0111axaxaxaxLnnnnnnnnnnnnnyyYaaAxxxxX001100,11在什么条件下）(0)det(X),1,0()(njyxLjjn)2(YXA 求ia三种插值三种插值方法方法有唯一解)2(1.1 1.1 拉格朗日插值多项式拉格朗日插值多项式nixxxxxxxxxxxxxxxxxlniiiiiiniii1,0,)()()()()()()(110110)3()()(0 xlyxLiniinjjnjiyxLjijixl)(,0,1)(则若又又（2）有唯一解，故有唯一解，故（3）与与（1）相同。相同。基函数基函数()ilx)1()(0111axaxaxaxLnn
6、nnn)2(YXA三种插值三种插值方法方法),(),()!1()()()()(0)1(baxxngxLxgxRnjjnnn1)1()(nnMg减小（粗略地看）如何使误差)(xRn平缓gjxx 接近njjnnxxnMxR01)!1()(三种插值三种插值方法方法1.2 1.2 误差估计误差估计增加n1.3 1.3 拉格朗日插值多项式的振荡拉格朗日插值多项式的振荡?)(?)(xRxLnnn55,11)(2xxxg63.363.3),()(limxxgxLnnRunge现象现象取n=2,4,6,8,10,计算Ln(x),画出图形-505-1.5-1-0.500.511.52y=1/(1+x2)n=2n
7、=4n=6n=8n=10三种插值三种插值方法方法Runge.m2.2.分段线性插值分段线性插值xjxj-1xj+1x0 xn其它,0,)()()(1111110jjjjjjjjjjjnjjjnxxxxxxxxxxxxxxxlxlyxI计算量与计算量与n n无关无关;n n越大，误差越小越大，误差越小.nnnxxxxgxI0),()(lim三种插值三种插值方法方法机翼下轮廓线3.3.三次样条插值三次样条插值样条函数的由来样条函数的由来飞机、船体、汽车外形等的放样（设计）飞机、船体、汽车外形等的放样（设计）细木条：样条细木条：样条3.3.三次样条插值三次样条插值,1,),()(1nixxxxsxS
8、iii,)()3),1,0()()2),1()()10223niiiiiiixxCxSniyxSnidxcxbxaxs数学样条（数学样条（spline）iiiidcban,4 个待定系数3）)1,1()()()()(),()(111 nixsxsxsxsxsxsiiiiiiiiiiii3）2），3）共 4n-2个方程三种插值三种插值方法方法自然边界条件）(0)()()40 nxSxS)(,)4)3)2xSdcbaiiii三次样条插值确定三次样条插值确定4n4n个系数需增加个系数需增加 2 2个条件个条件思考1）自然边界条件的几何意义是什么？2）样条插值为什么普遍用3次多项式，而不是2或4次？三
9、次样条三次样条插值插值).()(limxgxSn三种插值方法小结三种插值方法小结 拉格朗日插值（高次多项式插值）：拉格朗日插值（高次多项式插值）：曲线光滑；误差估计有表达式；收敛性不能保证。曲线光滑；误差估计有表达式；收敛性不能保证。用于理论分析，实际意义不大用于理论分析，实际意义不大。分段线性和三次样条插值（低次多项式插值）：分段线性和三次样条插值（低次多项式插值）：曲线不光滑（三次样条插值已大有改进）；误差估曲线不光滑（三次样条插值已大有改进）；误差估计较难（对三次样条插值）；收敛性有保证。计较难（对三次样条插值）；收敛性有保证。简单实用，应用广泛简单实用，应用广泛。其他：其他：Hermi
10、te插值、分段三次插值、二维插值等插值、分段三次插值、二维插值等根据需要，各取所需根据需要，各取所需。1.1.拉格朗日插值拉格朗日插值:自编程序自编程序,如名为如名为 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=inte
11、rp1(x0,y0,x,linear)3.3.三次样条插值三次样条插值:已有程序已有程序 y=interp1(x0,y0,x,spline)或或 y=spline(x0,y0,x)用MATLAB作插值计算注：注：MATLAB有样条工具箱（有样条工具箱（Spline Toolbox）用MATLAB作插值计算55,11)(2xxxg为例，作三种插值的比较为例，作三种插值的比较以以 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.3
12、077 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=
13、21个插值点，个插值点，三种方法作三种方法作插值，画图。插值，画图。chazhi1插值的应用加工时需要加工时需要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.051.0,2.50 1.2,2.051.4,1.69 1.6,1.40 1.8,1.182.0,1.00 2.2,0.862.4,0.74 2.6,0.64 模型模型 将图1逆时针方向转90度，轮廓线上下对称，只需对上半部计算一个
14、函数在插值点的值。图2 逆时针方向转90度的结果-5-4-3-2-101234500.511.522.533.544.55uv令v=x,u=-y 为什么要作数值积分为什么要作数值积分 许多函数许多函数“积不出来积不出来”,只能用数值方法，如只能用数值方法，如dxxxdxebabaxsin,22 积分是重要的数学工具，是微分方程、概率积分是重要的数学工具，是微分方程、概率论等的基础；在实际问题中有直接应用。论等的基础；在实际问题中有直接应用。对于用离散数据或者图形表示的函数对于用离散数据或者图形表示的函数，计算积分只有求助于数值方法。计算积分只有求助于数值方法。数值数值积分积分nabfIIdxx
15、fInkknnnba)(,lim)(1数数 值值 积积 分分 的的 基基 本本 思思 路路回回 忆忆 定定 积积 分分 的的 定定 义义各种数值积分方法研究的是各种数值积分方法研究的是k),(ba如何取值，区间如何取值，区间如何划分，如何划分，使得既能保证一定精度，计算量又小。使得既能保证一定精度，计算量又小。n n充分大时充分大时I In n就是就是I I的数值积分的数值积分（计算功效：算得准，算得快）（计算功效：算得准，算得快）1.1.从矩形公式到梯形公式数值积分数值积分yy=f(x)xbao)1(10nkknfhL)(,10kknkxffnabhbxxxxa)2(1nkknfhRnnRL
16、,平均，得到梯形公式)3()(2011nnkknffhfhTxk+1xkxk-1fk2.2.辛普森辛普森(Simpson)(Simpson)公式公式（抛物线公式）（抛物线公式）梯形公式相当于用分段线性插值函数分段线性插值函数代替)(xf每段要用相邻两小区间两小区间端点的三个函数值端点的三个函数值抛物线抛物线公式公式提高精度提高精度分段二次插值函数分段二次插值函数2221212222(,),(,),(,)0,1,1kkkkkkxfxfxfkm数值积分数值积分yy=f(x)xbaox2kf2kx2k+1x2k+2f2k+1f2k+2区间数必须为偶数区间数必须为偶数mn2)4(2),24(31121
17、01220mabhffffhSmkkmkkmm 对对k求和求和(共共m段段)，得（复合），得（复合）辛普森公式辛普森公式：)4(3)(22122222kkkxxkfffhdxxskk二次插值函数sk(x)构造用),(),(),(2222121222kkkkkkfxfxfx2.2.辛普森辛普森(Simpson)(Simpson)公式（抛物线公式）公式（抛物线公式）bannnTdxxfTITfR)(),(梯形公式在每小段上是用梯形公式在每小段上是用线性插值函数线性插值函数T T(x)代替代替 f(x),(,),)(2)()()(11 kkkkkkxxxxxxxfxTxf梯形公式梯形公式的误差估计的
18、误差估计)(2011nnkknffhfhTbadxxf)()(12)(2)()()(3111kxxkkkxxfhdxxxxxfdxxTxfkkkk 因为：因为：(x-xk)(x-xk+1)在在(xk,xk+1)不变号，所以：不变号，所以：)5()(12|),(|22abMhTfRn梯形公式梯形公式Tn的的误差是误差是h2阶的阶的),(,)(max2baxxfM 估计估计habn因为 103)(12|),(|nkknfhTfR 梯形公式梯形公式的误差的误差)()(121)(121 2afbfdxxfhTIban 103)(12)(nkkbannfhTdxxfTI)5()()(122afbfhTI
19、n同理可得：同理可得：)6()(180|),(|44abMhSfRn其中其中),(,)(max)4(4baxxfM辛普森公式辛普森公式Sn的误差是的误差是h4阶的阶的。辛普森公式的误差估计辛普森公式的误差估计梯形公式和辛普森公式的收敛性若对若对I某个数值积分某个数值积分In有有chIIpnnlim（非零常数）（非零常数）则称则称 In是是 p 阶收敛的阶收敛的。梯形公式梯形公式 2 2 阶收敛，辛普森公式阶收敛，辛普森公式 4 4 阶收敛。阶收敛。c=0:至少至少p阶收敛（超阶收敛（超p阶收敛）阶收敛）积分步长的自动选取积分步长的自动选取选定数值积分公式后，如何确定步长选定数值积分公式后，如何
20、确定步长h以满足给定的误差以满足给定的误差)()(122afbfhTIn梯形公式)(412nnTITInnTT2用二分法只要用二分法只要其中其中fk+1/2是原是原分点分点xk,xk+1的中点的中点（记记xk+1/2）的函数值的函数值1021222nkknnfhTT且且T T2n2n可在可在T Tn n基础上计算基础上计算)(3122nnnTTTInTI2)2(2nnhh高斯高斯(Gauss)(Gauss)求积公式求积公式矩形公式矩形公式(1)、(2)梯形公式梯形公式(3)辛普森公式辛普森公式(4)A Ak k是与是与f f无关的常数无关的常数代数代数精度精度设设,)(kxxf用用(7)计算计
21、算,)(badxxfI若对于若对于mk,1,0都有都有,IIn而当而当,1IImkn则则称称In的代数精度为的代数精度为m.)7()(1nkkknxfAINewton-Cotes方法方法梯形公式的代数精度（考察梯形公式的代数精度（考察T1）k=1f(x)=x222abxdxIba2)()()(21baabbfafhT3332abdxxIba2)(221baabTk=2f(x)=x2IT1IT 1梯形公式的代数精度为梯形公式的代数精度为1辛普森公式的代数精度为辛普森公式的代数精度为3高斯公式的思路高斯公式的思路取消对节点的限制，按照代数精度最大取消对节点的限制，按照代数精度最大的原则，同时确定节
22、点的原则，同时确定节点xk和系数和系数Ak构造求积公式构造求积公式)()(22112xfAxfAG对于对于11)(dxxfI使使G G2 2的代数精度为的代数精度为3 332,1)(xxxxf)()()(221111xfAxfAdxxf确定确定2121,AAxx03/202322311222211221121xAxAxAxAxAxAAA将将f(x)f(x)代入计算得代入计算得1,3/1,3/12121AAxx)3/1()3/1(2ffG用用n个节点，个节点，Gn的代数精度可达的代数精度可达2n-1,但是需解但是需解复杂的非线性方程组，实用价值不大。复杂的非线性方程组，实用价值不大。常 用 的
23、高 斯 公 式将将(a,ba,b)分小，把小区间变换为分小，把小区间变换为(-1(-1，1),1),再用再用G G2 2mkkkbazfzfhdxxf1)2()1()()(2)(322,3221)2(1)1(hxxzhxxzkkkkkkmkkhaxmabhk,1,0,/)(代数精度为代数精度为3节点加密时，原计算信息无法利用节点加密时，原计算信息无法利用思路思路：将积分区间分小，在小区间上用：将积分区间分小，在小区间上用n不太不太 大大的的 。而在节点加密一倍时能够利用原节点的函。而在节点加密一倍时能够利用原节点的函数值，可以把区间的端点作为固定节点。数值，可以把区间的端点作为固定节点。改进的
24、高斯公式nG)()()(121bfAxfAafAGnnkkknGauss-Lobatto求积公式求积公式 其中a,b为小区间的端点，nnAAxx,112为2n-2个参数，代数精度可达到代数精度可达到2n-3注意：实际计算中一般采用自适应方法确定步长注意：实际计算中一般采用自适应方法确定步长用用MATLAB 作数值积分作数值积分10nkknfhLnkknfhR1矩形矩形公式公式Sum(x)输入数组x(即fk),输出x的和(数)cumsum(x)输入数组x,输出x的依次累加和(数组)梯形梯形公式公式)(2011nnkknffhfhTtrapz(x)输入数组x,输出按梯形公式x的积分(单位步长)tr
25、apz(x,y)输入同长度数组 x,y,输出按梯形公式y对x的积分(步长不一定相等)用用MATLAB 作数值积分作数值积分mabhffffhSmkkmkkmn2),24(3112101220辛普森公式辛普森公式quad(fun,a,b,tol,trace)I,fn=quad()用自适应辛普森公式计算tol为绝对误差，缺省时为10-6Gauss-Lobatto公式公式)()()(121bfAxfAafAGnnkkknquadl(fun,a,b,tol,trace)I,fn=quadl()用自适应Gauss-Lobatto公式计算 tol为绝对误差，缺省时为10-6注意：注意：fun.m中应以自变
26、量为矩阵的形式输入中应以自变量为矩阵的形式输入(点运算点运算)矩形域上计算二重积分的命令：矩形域上计算二重积分的命令：dblquad(fun,xmin,xmax,ymin,ymax,tol)广义积分、二重和三重积分长方体上计算三重积分的命令：长方体上计算三重积分的命令：triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol)注：注：fun是被积函数，本身可以有自己的参数是被积函数，本身可以有自己的参数广义积分：广义积分：通过分析和控制误差，转换成普通积分通过分析和控制误差，转换成普通积分quadv(fun,a,b,tol,trace)向量值积分：向量值
27、积分：用用MATLAB 作数值积分作数值积分例例.计算计算4011s i nd xx1 1）矩形公式和梯形公式）矩形公式和梯形公式将将(0,(0,/4)100/4)100等分等分2 2）辛普森公式和辛普森公式和Gauss-Lobatto公式公式精确、方便精确、方便无法计算用数值给出的函数的积分无法计算用数值给出的函数的积分Jifen1a.mJifen1a.mJifen1b.mJifen1b.m精确值为精确值为2数值积分的应用数值积分的应用实例实例人造卫星轨道长度人造卫星轨道长度)20(sin,costtbytax决定由短半轴长半轴rssba,21dttbtadtyxL2022222022cos
28、sin44轨道长度轨道长度yxo 近地点s1=439km,远地点s2=2384kms1s2地球半径r=6371kmr需要作数值积分需要作数值积分s1=439km,s2=2384km,r=6371kmyxos1s2rs1s2yxoracb决定由短半轴长半轴rssba,21数值积分实例数值积分实例 人造卫星轨道长度人造卫星轨道长度dttbtaL202222cossin41222ssra7782.5212ssra1srac焦距212ssc7721.522cabdttbtaL202222cossin4用梯形公式和辛普森公式计算用梯形公式和辛普森公式计算只将区间只将区间5等分，梯形公式就给出很好的结果等
29、分，梯形公式就给出很好的结果轨道长度轨道长度 L L=4.8707 104千米千米数值积分实例数值积分实例 人造卫星轨道长度人造卫星轨道长度Jifen2.mJifen2.m布置实验布置实验目的目的1、掌握用、掌握用MATLAB计算拉格朗日、分段线性、三次计算拉格朗日、分段线性、三次样条三样条三 种插值的方法，改变节点的数目，对三种种插值的方法，改变节点的数目，对三种插值结果进行初步分析。插值结果进行初步分析。2、掌握用掌握用MATLAB及梯形公式、辛普森公式计算及梯形公式、辛普森公式计算数值积分。数值积分。3、通过实例学习用插值和数值积分解决实际问题。通过实例学习用插值和数值积分解决实际问题。内容内容实验练习实验练习 1010；1212；
