数值分析第8次课

1、1问题问题: 给定给定n个数据点个数据点 (xi , yi ) (i=1, 2, , n)ix2x1xnxiy1y2yny求直线求直线 y=a0+a1x 使得使得 niiixaay1102)(达到最小达到最小. 最小二乘一次多项式拟合最小二乘一次多项式拟合1 数据拟合的最小二乘法数据拟合的最小二乘法2 , , , , 21201011112xxxyxxyxxxxxyyxxxyyxxxxniiyiniixyniixniixxSnSSSnSaSnSSSSSaSSaaSSSnySyxSxSxS 设设 由上式求得由上式求得a0, a1, 代入代入 y=a0+a1x 得到得到最小二乘拟最小二乘拟合合(直

2、线直线)一次多项式一次多项式.3问题问题: 给定给定n个数据点个数据点 (xi , yi ) (i=1, 2, , n)ix2x1xnxiy1y2yny求求,2210 xaxaay 使得使得 niiiixaxaay122102)(达到最小达到最小. 最小二乘二次多项式拟合最小二乘二次多项式拟合4 niiiniiiniiniiniiniiniiniiniiniiniiyxyxyaaaxxxxxxxxn121121014131213121121 用用 Cholesky分解法求此对称正定阵分解法求此对称正定阵 用用 MATLAB 函数函数 z = Ar 由上式求得由上式求得a0, a1, a2, 得

3、到得到最小二乘拟合二次多项式最小二乘拟合二次多项式正则方程组正则方程组5xiyi例例 给定一组实验数据如下给定一组实验数据如下求求x, y的函数关系的函数关系.1 2 3 4 6 7 82 3 6 7 5 3 2u (1) 作散点分布图作散点分布图点的分布近似为抛物线点的分布近似为抛物线6 (2)确定近似表达式确定近似表达式设拟合曲线为二次多项式设拟合曲线为二次多项式22102)(xaxaaxPy (3) 建立正则方程组建立正则方程组, 7 n,3171 iix,179712 iix,1171713 iix,8147714 iix,2871 iiy,12171 iiiyx,635712 iii

4、xy7故正则方程组为故正则方程组为 6351212881471171179117117931179317210aaa (4) 求解正则方程组得求解正则方程组得,3864. 0,4318. 3,3182. 1310 aaa故所求拟合曲线为故所求拟合曲线为.3864. 04318. 33182. 1)(22xxxPy 8xiyi例例 给定一组实验数据如下给定一组实验数据如下求求x, y的函数关系的函数关系.1 2 3 4 6 7 82 3 6 7 5 3 2 Matlab解法解法: polyfit(1, 2, 3, 4, 6, 7, 8, 2, 3, 6, 7, 5, 3, 2, 2) ans=

5、-0.3864 3.4318 -1.31829例例 测得一发射源的发射强度测得一发射源的发射强度 I 与时间与时间 t 的一组数据如下的一组数据如下tiIi0.2 0.3 0.4 0.5 0.6 0.7 0.83.16 2.38 1.75 1.34 1.00 0.74 0.56试用最小二乘法确定试用最小二乘法确定 I 与与 t 的函数关系的函数关系.tIn (1)作散点分布图作散点分布图可以考虑用指数函数近似可以考虑用指数函数近似10 列数据表列数据表tiIi0.2 0.3 0.4 0.5 0.6 0.7 0.83.16 2.38 1.75 1.34 1.00 0.74 0.56lnIi1.1

6、506 0.8671 0.5596 0.2927 0 0.3011 0.5798求求lnI与与t的最小二乘直线的最小二乘直线. 将上表数据代入正则方程组将上表数据代入正则方程组得得 1858. 09891. 103. 25 . 35 . 3710aa其解为其解为.89. 2,73. 110 aa故所求拟合曲线为故所求拟合曲线为.64. 589. 289. 273. 1lnttIeeeI 11 列数据表列数据表tiIi0.2 0.3 0.4 0.5 0.6 0.7 0.83.16 2.38 1.75 1.34 1.00 0.74 0.56lnIi1.1506 0.8671 0.5596 0.29

7、27 0 0.3011 0.5798求求lnI与与t的最小二乘直线的最小二乘直线. 将上表数据代入正则方程组将上表数据代入正则方程组得得 1858. 09891. 103. 25 . 35 . 3710aa其解为其解为.89. 2,73. 110 aa故所求拟合曲线为故所求拟合曲线为.64. 589. 289. 273. 1lnttIeeeI Matlab解法解法:polyfit(0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.1506, 0.8671, 0.5596, 0.2927, 0, -0.3011, -0.5798, 1) ans= -2.8883 1.72

8、8312 求数据组的最小二乘拟合函数的步骤求数据组的最小二乘拟合函数的步骤 (1) 由给定数据确定近似函数的表达式由给定数据确定近似函数的表达式, 一般可一般可通过描点观察或经验估计得到通过描点观察或经验估计得到. (2) 按最小二乘原则确定表达式中的参数按最小二乘原则确定表达式中的参数, 即由即由残差平方和最小导出正则方程组残差平方和最小导出正则方程组, 求解得参数求解得参数.13第七章第七章 数值微分与数值积分数值微分与数值积分1 数值微分数值微分2 NewtonCotes求积公式求积公式3 复化求积公式复化求积公式4 Romberg求积公式求积公式14 利用离散点上函数的信息求函数导数近

9、似值利用离散点上函数的信息求函数导数近似值的方法的方法, 称为称为数值微分数值微分. 差商型数值微分公式差商型数值微分公式1 数值微分数值微分15由导数定义由导数定义hxfhxfxfh)()(lim)(0 当当h很小时很小时, 可用可用差商差商近似导数近似导数.16 差商型求导公式差商型求导公式 (3)中心差商公式中心差商公式0,)()()( hhxfhxfxf,)()()(hhxfxfxf .2)()()(hhxfhxfxf (1) 向前差商公式向前差商公式 (2) 向后差商公式向后差商公式17 几何意义几何意义hx xhx ABChxfhxfkBC)()( hhxfxfkAB)()( hh

10、xfhxfkAC2)()( B点切线斜率点切线斜率)(xf 从几何直观看从几何直观看: 中心差商效果最好中心差商效果最好18 截断误差截断误差)(2)( )()()( 1hOhhxfhxfhxfxf )(2)( )()()( 2hOhhxfhhxfxfxf )(12)()(2)()()( 223)3(3)3(hOhhxfhxfhhxfhxfxf 其中其中1,0321 由由Taylor公式可得公式可得19 二阶导数的中心差商公式二阶导数的中心差商公式2)()(2)()( hhxfxfhxfxf 截断误差截断误差)(12)()(2)()( )4(22 fhhhxfxfhxfxf 20近似计算近似计

11、算 badxxfI)(数值积分数值积分21 依据微积分基本定理依据微积分基本定理, 只要找到被积函数只要找到被积函数 f (x)的的原函数原函数 F (x), F (x)=f (x), 便有便有)()()(aFbFdxxfba 为什么还要对积分进行近似计算为什么还要对积分进行近似计算 大量的被积函数找不到用初等函数表示的原函数大量的被积函数找不到用初等函数表示的原函数 实验测量或数值计算给出的通常是一张离散函数实验测量或数值计算给出的通常是一张离散函数表表, 即被积函数的表达式未知即被积函数的表达式未知.数值积分数值积分222 Newton-Cotes 公式公式基本思想基本思想: 利用利用插值

12、多项式插值多项式).()(xfxLn banbadxxLdxxffI.)()()(其中其中Ln(x)是是n阶阶Lagrange插值多项式，用插值多项式，用Ln (x)的的积分近似积分近似 f (x)的积分，即的积分，即插值型求积公式插值型求积公式23 bankjjjkjkdxxxxxA,0)()(由由 决定决定,与与 无关无关.节点节点 f (x) 在在a, b上取上取 a x0 x1 0, 使得使得|,|max0knk 则称该求积公式是则称该求积公式是稳定稳定的的. 求积公式的稳定性求积公式的稳定性48 若求积公式是稳定的若求积公式是稳定的, 则则 f (x)的观察值的较小的的观察值的较小的

13、误差引起的求积结果的误差也是较小的误差引起的求积结果的误差也是较小的. 求积公式求积公式没有把没有把 f (x)的误差的误差“放大放大”很多很多.49), 1 , 0()(nkfxfkk nkkkknnfxfAfIfI0)()()(证明证明因此复化梯形公式是数值稳定的因此复化梯形公式是数值稳定的.当当).(2)1(2abhhnh 定理定理 复化梯形公式是复化梯形公式是数值稳定数值稳定的的. nkkkkfxfA0)( nkkA0 50 x0 x2xf (x)x4hhxn 2hxnmnnabh2, .hx3x1xn 1 复化复化Simpson公式公式分片二次多项式近似分片二次多项式近似51 将积分

14、区间将积分区间a, b划分为划分为n=2m等分等分, 步长步长 h=( b a )/n, 分点分点 xk= a+kh ( k=0, 1, , n). 在每个在每个小区间小区间 x2k 2 , x 2k ( k=1, , m)上用上用Simpson公式：公式： kkxxkkkkkxfxfxfxxdxxf222)()(4)(6)(21222222 )()(4)(321222kkkxfxfxfh 复化复化Simpson公式公式k=1, , m52 111122)(4)(2)()(3mkmkkkxfxfbfafh= Sn( f ) mkxxbadxxfdxxffIkk1222)()()( )()(4)

15、(3212221kkkmkxfxfxfh 复化复化Simpson公式公式53)()()(fSfIfRnn 当当 f (x)在在a, b上具有四阶连续导数时上具有四阶连续导数时, ),(),()(1)4()4(bamffmkk ),(222kkkxx 故得故得),(),(180)()(90)()4(4)4(5bafhabfmhfRn 复化复化Simpson公式的截断误差公式的截断误差 mkkfh1)4(5),(90 54 由复化由复化Simpson公式的截断误差知公式的截断误差知, 误差阶为误差阶为 h4, 收敛性是显然的收敛性是显然的, 事实上事实上, 只要只要 f (x) Ca, b则则可得

16、到可得到收敛性收敛性, 即即.)()(limdxxffSbann 由于求积系数均为正由于求积系数均为正, 与复化梯形公式一样的与复化梯形公式一样的证法可得复化证法可得复化 Simpson公式是公式是数值稳定数值稳定的的.55例例 对于函数对于函数,sin)(xxxf 给出给出n=8的函数表的函数表, 试用试用复化梯形公式及复化复化梯形公式及复化Simpson公式计算积分公式计算积分 10.sindxxxI解解ix)(ixf0.00.1250.250.3750.50.6250.750.8751.01.00.99739780.98961580.97672670.95885100.93615560.

17、90885160.84147090.8771925应用复化梯形公式求得应用复化梯形公式求得T8=0.9456909应用复化应用复化Simpson公式求得公式求得S8=0.9460832准确值准确值 I=0.9460831两者运算量基本相同两者运算量基本相同56trapz: 复化梯形公式求积分复化梯形公式求积分.用法用法: trapz(X, Y), 其中其中X, Y为相同维数的向量为相同维数的向量.例例: X=0.125:0.125:1.0;Y=sin(X)./X;X=0,X;Y=1,Y;trapz(X,Y)ans = 0.94569086358270Matlab函数函数57例例 若用复化求积公

18、式计算积分若用复化求积公式计算积分dxeIx 10的近似值的近似值, 若要求计算结果有若要求计算结果有4位有效数字位有效数字, n应取多大应取多大?解解, 1110 dxeIex.105 . 04 1 , 0, 1| )(|)( xexfxk复化梯形公式的误差复化梯形公式的误差)( )(12|2 fabhRT .83.40 n若用复化梯形公式求积分若用复化梯形公式求积分, n取取41能达到精度要求能达到精度要求.2121n 41021 58故应取故应取n=4. 该例表明该例表明, 为达到相同的精度为达到相同的精度, 用复化用复化Simpson公式所需的计算量比复化梯形公式要少公式所需的计算量比

19、复化梯形公式要少, 这也说明这也说明了复化了复化Simpson公式的精度高公式的精度高.复化复化Simpson公式的误差公式的误差)()(180|)4(4 fabhRS .25. 3 n41801n 41021 59 复化复化梯形公式的逐次分半算法梯形公式的逐次分半算法将区间将区间a, b分成分成n=2m等分等分, 记记,2mmabh .)(2)()(21212 mkmmkhafbfafhTm), 2 , 1 , 0( m称称 为为梯形值序列梯形值序列.2mT60 mhTTmm2122所有新增加节点的函数值之和所有新增加节点的函数值之和其中其中.2mmabh 复化复化梯形公式的逐次分半算法梯形

20、公式的逐次分半算法61以以n=8, m=3为例为例. 记记 fk= f (xk)x0 x2x4x6x3x1x5x7x8 )(2167654321808fffffffffabT )(21664280fffffab 75318ffffab 24T 所有新增加节点的函数值之和所有新增加节点的函数值之和. 3h62 复化梯形公式余项的后验估计复化梯形公式余项的后验估计);,(),(12112bafhabTIn );,(),(2122222bafhabTIn f ( 1 ), f ( 2 ) 分别是分别是 f (x) 在在a, b上的上的n个点与个点与 2n 个点处的算术平均值个点处的算术平均值 (每个

21、小区间上取一个点每个小区间上取一个点). 当当n较大时较大时, 有有.)( 1)()(21dxxfabffba 63因此因此, 若事先给定误差限若事先给定误差限 , 则当则当.3|2 nnTT时时, 就可停止计算就可停止计算, 并认为并认为 T2n是满足精度要求的近是满足精度要求的近似值似值.;412 nnTITI);()(21 ff ).(3122nnnTTTI 64 复化复化Simpson公式的逐次分半算法公式的逐次分半算法将区间将区间a, b分成分成 n=2m 等分等分, 记记,2mmabh evenodd2)(2)(4)()(3kmkmmkhafkhafbfafhSm, 2 , 1 m

22、称称 为为Simpson序列序列.2mS65;1612 nnSISI);,(),(18011)4(4bafhabSIn );,(),(218022)4(42bafhabSIn );()(2)4(1)4( ff ).(15122nnnSSSI 因此因此, 若事先给定误差限若事先给定误差限 , 则当则当.15|2 nnSS时可停止计算时可停止计算, 取取 S2n为满足精度要求的近似值为满足精度要求的近似值. 复化复化Simpson公式余项的后验估计公式余项的后验估计66nnnnnTTTTTI3134)(31222 4 Romberg求积公式求积公式启示启示: 是否用是否用 复化梯形公式余项的后验估

23、计表明复化梯形公式余项的后验估计表明nnTT31342 逼近逼近 I ( f ) 比用比用 T2n要好要好. 事实上有事实上有.313422nnnSTT 即梯形值序列的巧妙线性组合得到即梯形值序列的巧妙线性组合得到Simpson序列序列!67以以n=4为例加以说明为例加以说明. 记记 fk= f (xk), )(2167654321808fffffffffabT )(28642804fffffabT 8abkaxk )(422847654321808fffffffffabT .3)4(488TTS )( 2)( 48642753180fffffffffab 484TT 68).(15122nn

24、nSSSI .1516)(151222nnnnnSSSSSI 逼近逼近 I ( f ) 比用比用 S2n要好要好.回答回答: 是的是的, 记记,15)16(22nnnSSC 15)16(2nnSS 则它恰为复化则它恰为复化Cotes公式公式; 且有如下误差估计式且有如下误差估计式).(| )(|62hOfICn 复化复化Simpson公式余项的后验估计表明公式余项的后验估计表明 问题问题: 是否用是否用69.631636422nnnRCC 类似地可以得到类似地可以得到2nR其中其中被称为被称为Romberg序列序列.).(| )(|82hOfIRn 截断误差截断误差:701T2T4T8T16T

25、32T2S4S8S16S32S4C8C16C32C8R16R32R3422nnnTTS 151622nnnSSC 636422nnnCCR 停机准则停机准则:梯形值序列梯形值序列Simpson序列序列Cotes序列序列Romberg序列序列 Romberg求积公式求积公式71例例 计算计算.sin10dxxxI 2)1()0(1ffT =0.9207355)5 . 0(212112fTT =0.9397933 )43()41(412124ffTT=0.9445135=0.9456909 )87()85()83()81(812148ffffTT解解先求梯形值序列先求梯形值序列72nnTnSnCn

26、R24180.92073550.93979330.94451350.94569090.94614590.95608690.94608330.94608300.94608310.94608313422nnnTTS 151622nnnSSC .636422nnnCCR 利用只有两三位有效数字利用只有两三位有效数字的的T1, ,T8 经过三次外推得经过三次外推得到到7位有效数字位有效数字. 可见加速的可见加速的效果十分显著效果十分显著.用用Romberg算法计算如下算法计算如下73 理论依据理论依据: 复化梯形公式的余项展开复化梯形公式的余项展开.记记),(hTTn 定理定理 设设,)(baCxf 则则 kkhhhIhT24221)( 其中系数其中系数 k ( k=0, 1, )是与是与 h 无关的常数无关的常数. T (h) 逼近逼近 I 的速度是的速度是 O ( h2 )阶阶.74,3)()2(4)(262411 kkhhhIhThThT 当区间当区间a, b 2n等分时等分时, 则有则有),2(2hTTn 在定理中以在定理中以 h/2 代替代替 h 得得,2164224221 k