大学数学实验 实验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=1f(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
大学酸学实 前用的高斯公式 A1x1+A2x2=0 将f()代入计算得A1x2+42x=2/3 将(ab分小,把小区间变换为(-1,1),再用 A f(x)dx f(=k)+f(-k2 A=A2 x→(2=f(-1/5)+f/5 (1)_某+x-1 xk txk 用n个节点,G的代数精度可达2n1,但是需解 复杂的非线性方程组,实用价值不大 h=(b-a)/m, xk=a+ kh, k=0, 1,.m 代数精度为3 改进的高斯公式 用 MATLAB T值积分 思路:将积分区间分小,在小区间上用n不太大 矩形 的Gn。而在节 公式 数值,需要把 端点作为固定节点。可以改 进高斯公式 Sum(x)输入数组x即输出x的和(数 cumsum(x)输入数组x,输出x的依次累加和(数组) Gauss-Lobatto求积公式 梯形 G,=A, f(a)+2Af(x)+A,f(b) 公式 tapz(x)输入数组x,输出按梯形公式r的积分(单位步长 4,…A为2n2个参数 代数精度可达到2n-3 rap(xyl输入同长度数组x输出按梯形公式 y对x的积分(步长不一定相等) (大学数学实验) 用 MATLAB教积分 匚辛普森公式 重和三重积分 5n=3(6+/2m+42k+1+2)/ 2k), hs b-a 矩形域上计算二重积分的命令 k=0 Iquad( fun, xmin, xmax ymin, ymax, tol) I quad(@fun. n.a., b, tol)用自适应辛普森公式计算 tol为绝对误差,缺省时为105 fun是被积函数 Gauss.- Lobatto公式 G,=A,/(a)+2A/(2)+A,/(b) 长方体上计算三重积分的命令 dl(@fun, a, b, to)用自适应(asoh0式计算 triplequad(@ fun, xmin, xmax, ymin, ymax, zmin, zmax, tol) tol为绝对误差,缺省时为106 fun是被积函数
6 0 2 / 3 0 2 3 2 2 3 1 1 2 2 2 2 1 1 1 1 2 2 1 2 + = + = + = + = A x A x A x A x A x A x A A 将f(x)代入计算得 x1 = −1/ 3 , x2 = 1/ 3 , A1 = A2 = 1 ( 1/ 3) (1/ 3) 2 G = f − + f 用n个节点，Gn的代数精度可达2n-1, 但是需解 复杂的非线性方程组，实用价值不大。 常用的高斯公式 将(a,b)分小，把小区间变换为(-1，1), 再用G2 ∫ ∑= ≅ + m k k k b a f z f z h f x dx 1 (1) (2) [ ( ) ( )] 2 ( ) 2 2 3 , 2 2 3 (1) 1 (2) x x 1 h z x x h z k k k k k k − + + = + = − − h = (b − a)/ m, xk = a + kh, k = 0,1,L m 代数精度为3 思路：将积分区间分小，在小区间上用n不太 大 的 。而在节点加密一倍时能够利用原节点的函 数值，需要把区间的端点作为固定节点。可以改 进高斯公式 改进的高斯公式 Gn ( ) ( ) ( ) 1 2 G A1 f a A f x An f b n k n = + ∑ k k + − = Gauss-Lobatto求积公式 其中a, b为小区间的端点， n A An x2,L, x −1, 1,L 为2n-2个参数， 代数精度可达到2n-3 用MATLAB 作数值积分 ∑ − = = 1 0 n k n k L h f ∑= = n k n k R h f 1 矩形 公式 Sum(x) 输入数组x(即fk),输出x的和(数) cumsum(x) 输入数组x,输出x的依次累加和(数组) 梯形 公式 ( ) 2 0 1 1 n n k n k f f h T = h∑ f + + − = trapz(x) 输入数组x,输出按梯形公式x的积分(单位步长) trapz(x,y) 输入同长度数组 x,y,输出按梯形公式 y对x的积分(步长不一定相等) 用MATLAB 作数值积分 m b a f f f f h h S m k k m k n m k 2 ( 4 2 ), 3 1 1 2 1 0 0 2 2 1 − = + + ∑ + ∑ = − = − = + 辛普森公式 quad(@fun,a,b,tol) 用自适应辛普森公式计算 tol为绝对误差，缺省时为10-6 Gauss-Lobatto公式 ( ) ( ) ( ) 1 2 1 G A f a A f x A f b n n k n = + ∑ k k + − = 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是被积函数
用 MATLAB下教值积分 值积的应用 例计到 人造卫星轨道长度 sIn x 1)矩形公式和梯形公式 近地点s1=43%km,远地点s2=2384km 地球半径=6371km 将(0,π/4)100等分 长半轴 cost,y= bsin t b~短半轴 2)辛普森公式和 Gauss-Lobatto公式 (0≤t≤2n) 由s1,s2,r决定 精确 算用数值 函数的积分 L=42v22+jdt=412va2sin21+b cos2tdt 需要作数值积分 教值寂分臭人卫凰轨长度 积分例人鹭卫孰长度 L=4 用梯形公式和辛普森公式计算 s1=439km,2=2384km,r=6371km 2a=2r+52+5a=r+52+=782 轨道长度L=48707×104千米 焦距c=a-r-s1c=2 只将区间5等分,梯形公式就给出很好的结果 (学静学实鉴 布置实验 1、掌提用 MATLAB计算拉格朗日、分段线性、三次样 的方法,改变节点的数目,对三种插值 行初步分析 2、掌提用 MATLAB及梯形公式、辛普森公式计算数值 3、通过实例学习用播值和数值积分解决实际问题 内容课上布置。或见网络类堂
7 用MATLAB 作数值积分 例. 计算 4 0 1 1 sin d x x π − ∫ 1）矩形公式和梯形公式 将(0, π /4)100等分 2）辛普森公式和Gauss-Lobatto公式 精确、方便 无法计算用数值给出的函数的积分 MATLAB 5.3.lnk Jifen1.m 数值积分的应用 实例 人造卫星轨道长度 (0 2 ) cos , sin ≤ ≤ π = = t x a t y b t 由 决定 短半轴 长半轴 s s r b a , , ~ ~ 1 2 L x y dt a t b tdt ∫ ∫ = + = + 2 0 2 2 2 2 2 0 2 2 4 4 sin cos π π & & 轨道长度 y x o 近地点s1=439km,远地点s2= 2384km s s 1 2 地球半径r=6371km r 需要作数值积分 s1=439km, s2= 2384km, r=6371km y x o s s 1 2 r s1 s2 y x o r a c b a ~ 长半轴, b ~ 短半轴,由s1,s2 ,r决定 数值积分实例 人造卫星轨道长度 L a t b t dt ∫ = + 2 0 2 2 2 2 4 sin cos π 2 2 2 1 a = r + s + s 7782.5 2 2 1 = + = + s s a r 1 焦距c = a − r − s 2 2 1 s s c − = 7721.5 2 2 b = a − c = L a t b t dt ∫ = + 2 0 2 2 2 2 4 sin cos π 用梯形公式和辛普森公式计算 只将区间5等分，梯形公式就给出很好的结果 轨道长度 L=4.8707×104千米 数值积分实例 人造卫星轨道长度 MATLAB 5.3.lnk Jifen2.m 布置实验 目的 1、掌握用MATLAB计算拉格朗日、分段线性、三次样 条三 种插值的方法，改变节点的数目，对三种插值 结果进行初步分析。 2、 掌握用MATLAB及梯形公式、辛普森公式计算数值 积分。 3、 通过实例学习用插值和数值积分解决实际问题。 内容 课上布置，或参见网络学堂