大学数学实验 实验3的基本内容 插值的基本原理 Experiments in Mathematics 三种插值方法:拉格朗日插值,分段线性 插值,三次样条插值。 实验3插值与数值积分 2.插值的 MATLAB实现及插值的应用。 3.数值积分的梯形公式、辛普森公式和高斯公式 半大歇教科系 4.数值积分的 MATLAB实现及数值积分的应用。 什么是插值?从查函数表说起 少的越本原還 插值问题的提法 标准正态分布函数袭02=-∫e 已知n+1个节点(x1,y)(=0,1,…n,其中x 互不相同,不妨设a=x<x<…<xn=b) 100.84130.84380.8461 求任一插值点x(≠x)处的插值y 1.10.86430.86650.8686 节点可视为由 12|0.88490.886908888 y=g(x)产生, 求φ(1114) g表达式复杂 01141086+08065×0.4 甚至无表达式 值 (学静学实鉴 学学实纷 蜡的 本取求解问题的基本路 S三种1拉格朗日 Lagrange)多项式插值 1.0插值多项式 构造一个(相对简单的)函数y=f(x),通过全部节点即 L, (x=ax"+a,x+.+a, x+ao (1) f(x)=y,(j=0 再用f(x)计算插值,即y=f(x) Ln(x)=yG=0,1…n) XA=Y(2) de(x)≠0(在什么条件下):(2)有唯一解 xo x x
1 大学数学实验 Experiments in Mathematics 实验3 插值与数值积分 清华大学数学科学系 实验3的基本内容 3.数值积分的梯形公式、辛普森公式和高斯公式。 1.插值的基本原理； 三种插值方法：拉格朗日插 值，分段线性 插值，三次样条插值。 2.插值的 MATLAB 实现及插值的应用。 4.数值积分的 MATLAB 实现及数值积分的应用。 查函数表 什么是插值？从查函数表说起 ∫−∞ − Φ = x t x e dt 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个节点 ( x , y ) ( j 0,1, n, j j = L 其中 j x 互不相同，不妨设 ), 0 1 a x x x b = < <L< n = 求任一插值点 ( ) * j x ≠ x 处的插值 . * y • • • • • 0 x 1 x n x 0 y 1 y 节点可视为由 y = g ( x)产生, g表达式复杂, 甚至无表达式 * x * y • • • • • 0 x 1 x n x 0 y 1 y 求解插值问题的基本思路 构造一个(相对简单的)函数 y = f (x), 通过全部节点,即 f ( x ) y ( j 0,1, n ) j = j = L 再用 f(x)计算插值，即 ( ). * * y = f x * x * y 插值的 基本原理 1.拉格朗日(Lagrange)多项式插值 1.0 插值多项式 ( ) (1) 1 0 1 L x a x a 1x a x a n n n n = n + + + + − − L ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = − − n n n n n n n n y y Y a a A x x x x X M M L L L 0 0 1 1 0 0 , , 1 1 Q det( X ) ≠ 0 (在什么条件下） L (x ) y ( j 0,1, n) n j = j = L XA = Y (2) 求 i a 三种插值 方法 ∴(2)有唯一解
三种值 L(x)=a,r"+a-x+.+ax+ao 1.2误差估计 三种插 方油11拉格朗日播值多项式XA=Y(2) Ln(x)=∑y(x)(3) B()=g(3-4(=(m+71(-xe(cb l(x)= x-x)…(x-xx-x4)…(x-x) (-x)=0,n 8(5)sMalE /R,(x) s Mn+L[Tl (x-x)(x-x2 1(x1)=10.1≠j Ln(xj)=yj基函数 如何使误差|(x)减小(粗略地看) 平缓n增加 又(2)有唯一解,故(3)与(1)相同 1.3拉格朗日插值多项式的振荡 三种2.分欺能性 想n个→L,(x)?=|R,(x)↓ 5≤x≤5 取n=2,4,6,.8,10,计 算L(x),画出图形 ,x15x≤x,计算量与n无关 n越大,误差越小 liml(x)=g(x),x≤x≤xn Runge现象 lim L(x)=g(x),-363≤x≤363 (学静学实鉴 (大学数学实验) 3.三次样条插值 3.三次禅杂值[教学样条( spline)]「种 样条函数的由来 S(x)={s(x),x∈[x,x],i=1…n D)s(x)=ax+b 2)S(x)=y(=01,…n)4n个待定系数 3)S(x)∈CIx0,xn a1;b1,C;,d 细木条:样条 S;(x1)=S1+1(x1),s(x1)=s1+1(x) s:(x1)=s+1(x)(=1n-1) 飞机、船体、汽车外形等的放样(设计) 2),3’)共4n-2个方程
2 1.1 拉格朗日插值多项式 i n x x x x x x x x x x x x x x x x l x i i i i i i n i i n i L L L L L , 0,1 ( ) ( )( ) ( ) ( ) ( )( ) ( ) ( ) 0 1 1 0 1 1 = − − − − − − − − = − + − + ( ) ( ) (3) 0 L x y l x i n i n ∑ i = = i j n j j L x y i j i j l x ∴ = ⎩ ⎨ ⎧ ≠ = = ( ) 0, 1, Q ( ) 又（2）有唯一解，故（3）与（1）相同。 基函数 ( ) i l x ( ) (1) 1 0 1 1 L x a x a x a x a n n n n = n + + + + − − L XA=Y (2) 三种插值 方法 ( ), ( , ) ( 1)! ( ) ( ) ( ) ( ) 0 ( 1) x x a b n g R x g x L x n j j n n n − ∈ + = − = ∏= + ξ ξ 1 ( 1) ( ) + + ≤ n n g ξ M 如何使误差 Rn ( x) 减小（粗略地看） x 接近 x j g 平缓 ∏= + − + ≤ n j j n n x x n M R x 0 1 ( 1)! ( ) 三种插值 方法 1.2 误差估计 n增加 1.3 拉格朗日插值多项式的振荡 n ↑ ⇒ L ( x) ? ⇒ R ( x) ↓ ? n n , 5 5 1 1 ( ) 2 − ≤ ≤ + = x x g x lim ( ) = ( ), −3.63≤ ≤ 3.63 →∞ L x g x x n n Runge现象 取n=2,4,6,8,10,计 算Ln(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.分段线性插值 • • • • • • xj xj-1 xj+1 x0 xn ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ ⎨ ⎧ ≤ ≤ − − ≤ ≤ − − = = + + + − − − = ∑ 0 , 其它 , , ( ) ( ) ( ) 1 1 1 1 1 1 0 j j j j j j j j j j j n j n j j x x x x x x x x x x x x x x l x I x y l x 计算量与n无关; n越大，误差越小. n n n I x = g x x ≤ x ≤ x →∞ 0 lim ( ) ( ), 三种插值 方法 机翼下轮廓 线 3. 三次样条插值 样条函数的由来 飞机、船体、汽车外形等的放样（设计） 细木条：样条 3. 三次样条插值 ( ) { ( ), [ , ], 1, } S x = si x x∈ xi−1 xi i = Ln 3) ( ) [ , ] 2) ( ) ( 0,1, ) 1) ( ) ( 1, ) 0 2 3 2 n i i i i i i i S x C x x S x y i n s x a x b x c x d i n ∈ = = = + + + = L L 数学样条（spline） i i i di a b c n , , , 4 个待定系数 3） ( ) ( ) ( 1, 1) ( ) ( ), ( ) ( ) 1 1 1 ′′ = ′′ = − = ′ = ′ + + + s x s x i n s x s x s x s x i i i i i i i i i i i i 3’） 2），3’）共 4n-2个方程 三种插值 方法
大学酸学实 三种辑值力油小繪 三次禅辑值确4个系歙哪增加2个套件 拉格朗日插值〔高次多项式插值) 4)S"(x0)=S"(xn)=0(自然边界条件) 曲线光滑;误差估计有表达式;收敛性 不能保证(振荡现急) 2)3)4)→a1,b;,c1,d1→S(x) 用子论分析,安尾文不大。 limS(x)=g(x) ·分段线性和三次样条插值(低次多项式插值): 思考1)自然边界条件的几何意义是什么? 曲线不光滑(三次样条插值已大有改进);误差估 计较难(对三次样条插值);收敛性有保证 2)样条插值为什么普遍用3次多项式, 单用,用广促, 而不是2或4次? 用 MATLAB作插值计算 用 MATLAB作插值计算 1.拉格朗日插值:自编程序如名为1agx,m的M文件, 以g(+,-5≤x≤5为例,作三种插值的比较 第一行为 function y=1agr(x0,y0,x) 010M1品 输入:节点x0,y0,插值点x(均为数组,长度自定义)) 用n=11个节 0.50000.80000.84340.75000.82 输出:播值y(与x同长度数组) 点,m=21个1000.50000.50000.50000.5000 应用时输入x0,y0,x后,运行y=1agx(x0,y 插值点,三 种方法作插 2.00000.2000 分段线性插值:已有程序y= interp1(x0,y0,x) 值,画图。 0.15000.1401 3.00000.1000 0.10000.1000 3.三次样条插值:已有程序y= interp1(x0,y0,x,'sp1ine") 或y=sp1ine(x0,yo,x) 4.00000.05880.05880.05880.0588 注:1agxm程序可参考课本 hazel 4.50000.04711.57870.04860.0484 5.00000.03850.0385003850.0385 (学酸学实) 插值的应用 数控机床加工零件 人值为什么要作死学x 表1间隔02的加如工坐标x(图1右半部的数据) ·积分是重要的数学工具,是微分方程、概率 D0s0024710443106368083051 1025012205141691614018118 论等的基础;在实际问题中有直接应用。 2010028414074260644m5的数 ·许多函数“积不出来”,只能用数值方法,如 模型将图1逆时针方向转90度,轮廓线上下 加工时警要x每对称,只需对上半部计算一个函数在插值点 改变0.05时的y值的值 对于用离散数据或者图形表示的函数, 计算积分只有求助于数值方法
3 4) S ′′( x0 ) = S ′′( xn ) = 0 (自然边界条件） 2 ) 3 ) 4 ) a , b , c , d S ( x ) ⇒ i i i i ⇒ 三次样条插值确定4n个系数需增加 2个条件 思考 1）自然边界条件的几何意义是什么？ 2）样条插值为什么普遍用3次多项式， 而不是2或4次？ 三次样条 插值 limS(x) g(x). 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作插值计算 , 5 5 1 1 ( ) 2 − ≤ ≤ + = x x 以 g x 为例，作三种插值的比较 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.0,1.00 2.2,0.86 2.4,0.74 2.6,0.64 ……… 1.0,2.50 1.2,2.05 1.4,1.69 1.6,1.40 1.8,1.18 0.0,5.00 0.2,4.71 0.4,4.31 0.6,3.68 0.8,3.05 模型 将图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 e dx b a b a x ∫ ∫ − sin , 2 2 • 积分是重要的数学工具，是微分方程、概率 论等的基础；在实际问题中有直接应用。 • 对于用离散数据或者图形表示的函数， 计算积分只有求助于数值方法。 数值 积分
大学静学实扮 款积分的基本路 少值积分1.从矩形公式到梯形公式 y=x)a=1<x1<…Ik<xn=b 回钇积分的宠文 fk= f(K) I=f()dx= lim In, In=>f(Ek) (1) 分大时J就是酌的数值积分 各种数值积分方法研究的是 (2) 5k如何取值,区间(a,b)如何划分, Ln,Rn平均,得到 使得既能保证一定精度,计算量又小 梯形公式Tn=b>+(6+1m)(3) 数值积分2.辛普森( Simpson)公式 少2.辛普森( Simpson)公式(抛物线公式) (抛物线公式) 用xx,x)(x1,4)、(x2x2,2x2)构造二次插值函数s(x) 梯形公式相当于用分段线性插值函数代替f(x) 提高精度]分段二次播值函数口魏物线 Sk(x)x=(2k+4/2k+1+/2k+2) 每段要用相邻两小区间 对k求和(共m段),得辛普森公式 端点的三个函数值 区间数必须为偶数n=2m a xx x+ IN+, b (xxf)(x,x)(x+2,f2) Sm=(+12m+4/2++2∑/众b-0(可 k=0,1…;m-1 学学实 形公T=b+;(+f)≈f(x)d f(x)d b∑r"(5,) 的惧计 I-T R(f,T。) f(x)d T 12」f(x) (f(b)-f'(a 梯形公式在每小段上是用线性插值函数x)代替f(x) 梯形公式 的误R(f,Tn)|≤ ∑|"(5 (x)=7(x)+/(0x-x)x-x1)xn∈(x,x1) 估计M2=maf'(x)x∈(a,b)因为m=a 因为:(xxk)(x-x在(xx4+1)不变号,所以: U()-T(xdr =/152 R(f,7)1,M2(b-a)(5)命梯形公式厂的 2](x-Xx+1)=-12f) 差是2阶的
4 n b a I f x dx I I f n k n n k n b a − = ∫ = =∑= →∞ ( ) lim , ( ) 1 ξ 数值积分的基本思路 回忆定积分的定义 各种数值积分方法研究的是 ξ k 如何取值，区间 如何划分， (a,b) 使得既能保证一定精度，计算量又小。 n充分大时In就是I的数值积分 数值积分 1.从矩形公式到梯形公式 y y=f(x) b x o a (1 ) 1 0 ∑ − = = n k n k L h f , ( ) , 0 1 k k k n f f x n b a h a x x x x b = − = = < < L L < = ( 2 ) 1 ∑= = n k n k R h f Ln R n , 平均，得到 梯形公式 ( ) (3) 2 0 1 1 n n k n k f f h T = h∑f + + − = xk+1 x xk-1 k fk 2.辛普森(Simpson)公式 （抛物线公式） 梯形公式相当于用分段线性插值函数代替 f (x) 每段要用相邻两小区间 端点的三个函数值 抛物线 公式 提高精度 分段二次插值函数 2 2 21 21 22 22 ( , ),( , ),( , ) 0,1, , 1 kk k k k k xf x f x f k m ++ + + = − L 数值积分 y y=f(x) b x o a x2k f2k x2k+1 x2k+2 f2k+1 f2k+2 区间数必须为偶数 n = 2m (4) 2 ( 4 2 ), 3 1 1 2 1 0 0 2 2 1 m b a f f f f h h S m k k m k m m k − = + + ∑ + ∑ = − = − = + 对k求和(共m段)，得辛普森公式： ( 4 ) 3 ( ) 2 2 1 2 2 2 2 2 = + + + + ∫ + k k k x x k f f f h s x dx k k 二次插值函数sk 用( , ),( , ),( , )构造 (x) 2k 2k 2k+1 2k+1 2k+2 2k+2 x f x f x f 2.辛普森(Simpson)公式（抛物线公式） ∫ = − b a n T n R ( f , T ) f ( x ) dx 梯形公式在每小段上是用线性插值函数T(x)代替 f(x) ( )( ), , ( , ) 2 ( ) ( ) ( ) − − +1 ∈ +1 ′′ = + k k k k k k x x x x x x x f f x T x η η 梯形公式 的误差估计 ( ) 2 0 1 1 n n k n k f f h T = h∑ f + + − = ∫ ≈ b a f (x)dx ( ) 12 ( )( ) 2 ( ) [ ( ) ( )] 3 1 1 1 k x x k k k x x f h x x x x dx f f x T x dx k k k k ξ ξ − − = − ′′ ′′ − = ∫ ∫ + + + 因为：(x-xk)(x-xk+1)在(xk,xk+1)不变号，所以： ( ) (5) 12 | ( , ) | 2 2 M b a h R f Tn ≤ − 即梯形公式Tn的 误差是h2阶的 max ( ) , ( , ) 2 估计 M = f ′′ x x∈ a b h b a n − 因为 = ∑ − = ≤ ′′ 1 0 3 ( ) 12 | ( , ) | n k n k f h R f T ξ 梯形公式 的误差 ( ' ( ) ' ( )) 12 1 ( ) 12 1 ( ) 12 ( ) ' ' 2 1 0 3 f x dx f b f a h I T f h I T f x dx T b a n n k k b a n n → − = − − − − = − = − ′′ ∫ ∫ ∑ − = ξ
大学酸学实 辛普森公式的误差估计 梯形公式和辛普森公式的收做性 同理可得 若对l某个数值积分有m n少0hPC(非零常数 R(f,S)M4(b-a)(6 则称l是p阶收敛的 180 →梯形公式2阶收斂,辛普森公式4阶收敢。 其中M4=ma|r(4(xxe(a,b) 即辛普森公式S的误差是h4阶的 积分步长的自动选取 高斯( Gauss)求积公式 选定敷值积分公式后,如何确定步长以满足给定的误差g 矩形公式(1)、(2 Akf(xk)(7) 梯形公式1-.-1(f(b)-f(a)b→%(→2n 梯形公式 1-7n=(1-T)口1-T2n≈3(72m-T 方法辛普森公式(4 是与关的常数 用二分法只要2n-则≤E曰|-2 设f(x)=x,用(计算 f(x)dx 且Tm可在Tn 基础上计算72n=2 n+∑f 若对于k=0,1,…m都有Jn=l, 其中f2是原分点xx+,的中点(记x+m2)的函数值 而当k=m+1,n≠1,则称J的代数精度为m (学静学实鉴 (大学数学实验) 梯形公式的代教精度(考察T1) 高斯公式的尾路 T=l(a)+f(b) b-a)(a+b) 取消对节点的限制,按照代教精度最大 的原则,同时确定节点x和系数A f(x)x 对于I f(x)dx构造求积公式 G2=A1f(x1)+A2f(x2)使G2的代数精度为3 f(x)=1x2,x3|fx)hk=4f(x1)+2f(x2) 梯形公式的代数精度为1辛普森公式的代数精度为3 确定x1,x2,A1,A2
5 同理可得： ( ) (6) 180 | ( , )| 4 4 M b a h R f Sn ≤ − 其中 max ( ) , ( , ) (4) M 4 = f x x ∈ a b 即辛普森公式Sn的误差是h4阶的。 辛普森公式的误差估计 梯形公式和辛普森公式的收敛性 若对I某个数值积分In有 c h I I p n n = − →∞ lim （非零常数） 则称 In是 p 阶收敛的。 梯形公式 2 阶收敛，辛普森公式 4 阶收敛。 积分步长的自动选取 选定数值积分公式后，如何确定步长h以满足给定的误差ε ( ( ) ( )) 12 ' ' 2 f b f a h 梯形公式 I − Tn ≈ − − ( ) 4 1 2 n Tn I − T ≈ I − − ≤ ε T2n Tn 用二分法只要 其中fk+1/2是原分点xk,xk+1的中点（记xk+1/2）的函数值 ∑ − = + = + 1 0 2 2 1 2 2 n k k n n f T h T 且T2n可在Tn 基础上计算 ( ) 3 1 T2 n T2 n Tn I − ≈ − − ≤ ε T n I 2 ( 2 ) 2 n n h h → → 高斯(Gauss)求积公式 矩形公式(1)、(2) 梯形公式(3) 辛普森公式(4) Ak是与f无关的常数 代数 精度 设 ( ) , k f x = x 用(7)计算 ( ) , ∫ = b a I f x dx 若对于 k = 0 ,1, L m 都有 I I, n = 而当 k m 1, I I, = + n ≠ 则称In的代数精度为m. ( ) (7) 1 ∑= = n k n k k I A f x Newton -Cotes 方法 梯形公式的代数精度（考察T1） k=1 f(x)=x 2 2 2 b a I xdx b a − = = ∫ 2 ( )( ) [ ( ) ( )] 2 1 b a a b f a f b h T − + = + = 3 3 3 2 b a I x dx b a − = = ∫ 2 ( )( ) 2 2 1 b a a b T − + = k=2 f(x)=x2 T = I 1 T ≠ I 1 梯形公式的代数精度为1 辛普森公式的代数精度为3 高斯公式的思路 取消对节点的限制，按照代数精度最大 的原则，同时确定节点xk和系数Ak 构造求积公式 ( ) ( ) 2 1 1 2 2 G = A f x + A f x 对于 ∫− = 1 1 I f ( x ) dx 使G2的代数精度为3 2 3 f (x) =1, x, x , x ( ) ( ) ( ) 1 1 2 2 1 1 f x dx = A f x + A f x ∫− 确定 1 2 1 2 x , x , A , A