《数值分析与算法》第七讲数值积分与数值微分

数值分析(7)NumericalAnalysisWenjianYu2第七章数值积分与数值微分WenjianYu3数值积分的基本概念WenjianYu4数值积分目的与用途经典问题:算几何形体的面积、体积，力学中物体的重心位置例:铝制波纹瓦的长度问题由一块平整的铝板压制而成.若每个波纹的高度(自中心线)为1英寸,周期为2英寸,做4英尺长波纹瓦需多长铝板?第二类椭圆积分,无法解析求出!WenjianYu5数值积分基本思想...()fx积分系数积分节点希望用较少的计算量得到较准确的结果WenjianYu6插值型求积公式中矩形公式梯形公式WenjianYu7积分余项与代数精度反映了计算的截断误差插值余项的积分衡量求积公式准确度的另一个指标注意:对某些情况，代数精度并不是越高越好WenjianYu8积分余项与代数精度(至少0次代数精度)它至少有n次代数精度即插值型可推广WenjianYu9积分余项与代数精度WenjianYu10求积公式的收敛性与稳定性(一系列求积公式的性质)积分问题一般不太敏感WenjianYu11求积公式的收敛性与稳定性这是控制数值计算误差能达到的最佳情况要尽量寻求稳定的公式(估计截断误差)WenjianYu12牛顿-柯特斯公式WenjianYu13Newton-Cotes公式这就是n阶牛顿-柯特斯公式n=1,1/2,1/2n=2,1/6,2/3,1/6n=4,7/90,16/45,2/15,16/45,7/90Cotes系数表•一系列求积公式•便于使用WenjianYu14Newton-Cotes公式n=1,1/2,1/2n=2,1/6,2/3,1/6n=4,7/90,16/45,2/15,16/45,7/90n=8,Cotes系数表思考题梯形公式Simpson公式Cotes公式中矩形公式可看成是n=0时的特例WenjianYu15Newton-Cotes公式(最准确)WenjianYu16Newton-Cotes公式关键看积分:(n阶公式至少有n次代数精度)一般不用n=3对应的N-C公式WenjianYu17低阶N-C公式的积分余项不保号,无法用积分中值定理2详细过程看书pp.246WenjianYu18稳定性、收敛性n=8,实际只使用n8的偶数阶N-C公式(代数精度不是越高越好)WenjianYu19复合求积公式WenjianYu20复合求积公式(compositequadrature)积分误差:n增大,误差减小仍是“机械求积公式”WenjianYu21复合求积公式2阶准确度WenjianYu22复合求积公式与复合梯形公式对比,看例7.4WenjianYu23复合求积公式步长折半的复合求积公式计算积分余项公式包含被积函数的高阶导数,很难应用.常常动态地确定步长h常用的动态减小步长策略是:步长折半,利用已算出的结果复合梯形公式的情况递推化的复合梯形公式:(逐渐减小,直到满足精度要求)只需再计算新增节点的函数值WenjianYu24复合求积公式步长折半的复合求积公式计算复合Simpson公式的情况很少使用中矩形公式的原因与梯形公式有相同的代数精度/准确度,计算量更小可类似构造复合中矩形公式,但在步长折半时,无法重用以前的结果WenjianYu25Remberg积分算法WenjianYu26复合梯形公式的余项展开式Th7.5所有小区间的积分求和:WenjianYu27复合梯形公式的余项展开式Th7.5所有小区间乘h/2求和:WenjianYu28RichardsonExtrapolation(“0”代表未经外推的原始公式)更准确的值!WenjianYu29RichardsonExtrapolationWenjianYu30Romberg算法龙贝格算法列三角形表格,按行依次计算计算公式可证明:具有2k+1次代数精度(类似高阶差商的计算)WenjianYu31Romberg积分算法基于等距分布的积分节点,得到稳定、准确的积分结果要求被积函数充分光滑!WenjianYu32Romberg算法011/80.99739781/40.98961583/80.97672671/20.95885115/80.93615563/40.90885177/80.877192610.841471010.92073551/20.93979330.94614591/220.94451350.94608690.946083001/230.94608330WenjianYu33Romberg算法10.5000001/20.4267770.4023691/220.4070180.4004320.4003021/230.4018120.4000770.4000540.4000501/240.4004630.4000140.4000090.4000090.4000091/250.4001180.4000020.4000020.4000020.4000020.400002准确值为2/5WenjianYu34自适应积分算法WenjianYu35自适应积分算法(adaptivequadrature)基本思想Romberg算法效果不好的情况积分节点没必要均匀分布怎样自动地非均匀取点,使计算结果达到精度要求?1.评估当前区间积分结果的准确度,若不准确就将它折半,直至小区间的结果准确2.用两个不同的求积公式算同一个积分,它们之差可近似判断结果的准确度(递归计算过程)Simpson公式与复合Simpson公式WenjianYu36自适应积分算法一个自适应求积算法对每个区间,用Simpson公式、复合Simpson公式计算无论区间大小,用相同的误差阈值;函数的递归调用原理算法:实际的算法需保证函数值不重复计算见课本的quadtx程序演示模块7.3,quadguiWenjianYu37自适应积分算法一个自适应求积程序quadtx例子:更多讨论通过阈值设置控制相对误差;不连续函数的特殊处理还有其他估计积分误差的方法,比如利用中矩形公式,梯形公式的差注意与Romberg算法的不同(取33个、65个点分别算出30.0548,29.8474,不如自适应积分准)102030405060708090100Q=29.8583,fcount=69[Q,fcnt]=quadtx(@humps,0,1,1e-3)fcnt=69,用了69个积分点WenjianYu38高斯求积公式WenjianYu39高斯求积公式解得:WenjianYu40高斯求积公式高斯积分有2n+1次代数精度插值型求积公式、代数精度的概念也可扩展(Th7.9)WenjianYu41高斯求积公式Th7.7比自适应积分算法使用方便pp.195性质(5)WenjianYu42高斯-勒让德公式012345高斯-勒让德积分表WenjianYu43高斯-勒让德公式WenjianYu44Matlab中的积分计算Matlab中的数值积分指定被积函数匿名函数@输入参数可以是向量，因此需采用逐项运算符号M文件:可处理含奇异点的积分含多个自变量:(第1个是主自变量，其他为参数)14011dxxf=@(x)1./sqrt(1+x.^4)0sinxdxxFunctionf=sinc(r)ifx==0f=1;elsef=sin(x)./x;endf_beta=@(t,z,w)t.^(z-1).*(1-t).^(w-1);1110(,)(1)zwzwttdt函数句柄@sincWenjianYu45Matlab中的数值积分一维积分的命令quad(第7.5节的自适应积分)quadl,quadgk(扩展的Gauss自适应积分,pp.263)Q=quad(f,0,1)quad(@sinc,0,pi)%functiondefinedbyM-fileans=0.58949[beta2_5,fcnt]=quad(f_beta,0,1,1e-5,0,2,5)beta2_5=0.033333fcnt=17[q,fcnt]=quad(fun,a,b,tol,trace,p1,p2,...)积分参数是否输出函数计算次数等.0/非0准确度控制：缺省值10-6WenjianYu46Matlab中的数值积分二重、三重积分dblquad,quad2d,triplequad符号积分定义符号变量:sym(),syms,符号积分:intsimple(表达式化简)double离散数据点积分复合梯形法trapz(x,y)symsxh=1/((x-.3)^2+.01)+1/((x-.9)^2+.04)-6h=1/((x-3/10)^2+1/100)+1/((x-9/10)^2+1/25)-6I=int(h)%不定积分I=10*atan(10*x-3)+5*atan(5*x-9/2)-6*xD=simple(int(h,0,1))%定积分D=5*atan(16/13)-6+10*piQexact=double(D)Qexact=29.858不用逐项运算符Matlab演示WenjianYu47WenjianYu48数值微分WenjianYu49数值微分(向前差分)(向后差分)(中心差分)利用Taylor展开推出:2阶准确度WenjianYu50数值微分(截断误差~舍入误差)此误差限何时最小?的误差限0.0754与实验相符!WenjianYu51数值微分WenjianYu52数值微分h0.10.45160490810.050.45407616940.45489992310.0250.45469262880.45489811520.4548979947p359准确的有效数字位数WenjianYu53数值微分的应用二阶中心差分例如n=2对每个格点列,微分方程变为代数方程