数值计算的资料

数值计算的资料第1页/共51页WenjianYu2第七章数值积分与数值微分

第2页/共51页3数值积分的基本概念第3页/共51页4数值积分目的与用途经典问题:算几何形体的面积、

体积，力学中

物体的重心位置例:铝制波纹瓦的长度问题由一块平整的铝板压制而成.

若每个波纹的高度(自中心线)

为1英寸,周期为2英寸,做4英尺长波纹瓦需多长铝板?

第二类椭圆积分,无法解析求出!

第4页/共51页5数值积分基本思想

...

积分系数

积分节点希望用较少的计算量得到较准确的结果第5页/共51页6插值型求积公式

中矩形公式梯形公式第6页/共51页7积分余项与代数精度

反映了计算的截断误差插值余项的积分衡量求积公式准确度的另一个指标注意:对某些情况，代数精度并不是越高越好第7页/共51页8积分余项与代数精度

(至少0次代数精度)它至少有n次代数精度即插值型

求积公式第8页/共51页9求积公式的收敛性与稳定性

(一系列求积公式的性质)

积分问题一般不太敏感

第9页/共51页10求积公式的收敛性与稳定性

这是控制数值计算误差能达到的最佳情况

要尽量寻求稳定的求积公式第10页/共51页11牛顿-柯特斯公式第11页/共51页12Newton-Cotes公式

就是n阶牛顿-柯特斯公式

n=1,1/2,1/2n=2,1/6,2/3,1/6

n=4,7/90,16/45,2/15,16/45,7/90Cotes系数表

一系列求积公式便于使用第12页/共51页13Newton-Cotes公式

n=1,1/2,1/2n=2,1/6,2/3,1/6

n=4,7/90,16/45,2/15,16/45,7/90

n=8,Cotes系数表思考题

梯形公式Simpson公式Cotes公式中矩形公式可看成是n=0时的特例.例题(板书)第13页/共51页14Newton-Cotes公式

关键看积分:

(n阶公式至少有n次代数精度)一般不用n=3对应的N-C公式第14页/共51页15低阶N-C公式的积分余项

不保号,无法用积分中值定理

2

详细过程自己看书第15页/共51页16稳定性、收敛性

n=8,实际只使用n<8的偶数阶N-C公式,也看出代数精度不是越高越好第16页/共51页17复合求积公式第17页/共51页18复合求积公式

(compositequadrature)

积分误差:n增大,误差减小

仍是机械求积公式第18页/共51页19复合求积公式

2阶准确度

第19页/共51页20复合求积公式

与复合梯形公式对比,例7.4第20页/共51页21复合求积公式步长折半的复合求积公式计算积分余项公式包含被积函数的高阶导数,很难应用.

常常动态地确定步长h常用的动态减小步长策略是:步长折半,利用已算出结果复合梯形公式的情况递推化的复合梯形公式:

(逐渐减小,直到满足精度要求)

只需再计算新增节点的函数值第21页/共51页22复合求积公式步长折半的复合求积公式计算复合Simpson公式的情况很少使用中矩形公式的原因

与梯形公式有相同的代数精度/准确度,计算量更小

可类似构造复合中矩形公式,但在步长折半时,无法重用以前的结果第22页/共51页23Remberg积分算法第23页/共51页24复合梯形公式的余项展开式

(可与Richardson外推结合)Th7.5

所有小区间的积分求和:第24页/共51页25复合梯形公式的余项展开式

Th7.5

所有小区间乘h/2求和:

第25页/共51页26RichardsonExtrapolation

(“0”代表未经外推的原始公式)

第26页/共51页27RichardsonExtrapolation

第27页/共51页28Romberg算法龙贝格算法列三角形表格,按行依次计算计算公式

可证明:

具有2k+1次代数精度

(类似高阶差商的计算)第28页/共51页29Romberg算法

做等距的节点分布,充分利用其上函数值得到最高准确度的结果要求被积函数充分光滑!第29页/共51页30Romberg算法

011/80.99739781/40.98961583/80.97672671/20.95885115/80.93615563/40.90885177/80.877192610.841471010.9207355

1/20.93979330.9461459

1/220.94451350.94608690.946083001/230.9460833第30页/共51页31Romberg算法

10.500000

1/20.4267770.402369

1/220.4070180.4004320.400302

1/230.4018120.4000770.4000540.400050

1/240.4004630.4000140.4000090.4000090.400009

1/250.4001180.4000020.4000020.4000020.4000020.400002

第31页/共51页32自适应积分算法第32页/共51页33自适应积分算法(adaptivequadrature)基本思想Romberg算法效果不好的情况积分节点没必要均匀分布怎样自动地非均匀取点,使计算

结果达到精度要求?1.评估当前区间积分结果的准

确度,若不准确就将它一分为二,直至小区间的结果准确2.

用两个低阶求积公式算同一个积分,它们之差可近似判断结果的准确度

(递归计算过程)第33页/共51页34自适应积分算法一个自适应求积算法对每个区间,用Simpson公式、复合Simpson公式计算无论区间大小,用相同的误差阈值;函数的递归调用原理算法:

实际的算法需保证函数值不重复计算见课本的quadtx程序

演示模块7.3,quadgui第34页/共51页35自适应积分算法一个自适应求积程序quadtx例子:更多讨论通过阈值设置控制相对误

差;不连续函数的特殊处理还有其他估计积分误差的

方法,比如利用中矩形公式,梯形公式的差注意与Romberg算法的不同

[Q,fcnt]=quadtx(@humps,0,1,1e-3)fcnt=69,用了69个积分点第35页/共51页36高斯求积公式第36页/共51页37高斯求积公式

解得:

第37页/共51页38高斯求积公式

高斯积分有2n+1次代数精度

插值型求积公式,代数精度的概念也可扩展

第38页/共51页39高斯求积公式

Th7.7比自适应积分算法使用方便第39页/共51页40高斯-勒让德公式

012345高斯-勒让德积分表第40页/共51页41高斯-勒让德公式

第41页/共51页42数值微分第42页/共51页43数值微分

(向前差分)(向后差分)(中心差分)利用Taylor展开推出:

第43页/共51页44数值微分

第44页/共51页45数值微分

h0.10.4516049081

0.050.45407616940.4548999231

00250.45469262880.45489811520.4548979947p359

准确的有效数字位数第45页/共51页46数值微分的应用

二阶中心差分

例如n=2

第46页/共51页47Matlab中的积分计算第47页/共51页Matlab中的数值积分指定被积函数inline命令Matlabv7推荐使用匿名函数@输入参数可以是向量，因此需采用逐项运算符号M文件:可处理含奇异点的积分带多个参数:>>f1=inline('1./sqrt(1+x.^4)'

)>>f2=@(x)1./sqrt(1+x.^4)Functionf=sinc(r)ifx==0f=1;elsef=sin(x)./x;end>>f_beta=inline('t.^(z-1).*(1-t).^(w-1)','t','z','w')函数句柄第48页/共51页Matlab中的数值积分一维积分的命令quad(第7.5节的自适应积分)quadl,quadgk(扩展的Gauss自适应积分,pp.231)>>Q=quad(f1,0,1)

%inlinefunction>>quad(@sinc,0,pi)%functiondefinedbyM-fileans=0.58949>>[beta2_5,fcnt]=quad(f_beta,0,1,1e-5,0,2,5)beta2_5=0.033333

fcnt=17[q,fcnt]=quad(fun,a,b,tol,trace,p1,p2,...)积分参数是否输出函数计算次数等.0/非0准确度控制：缺省值10-6第49页/共51页Matlab中的数值积分二重、三重积分dblquad,quad2d,triplequad符号积分定义符号变量:

sym(),syms,符号积分:intsimple(表达式化简)double离散数据点积分复合梯形法trapz(x,y)>>symsx>>h=1/((x-.3)^2+.01)+1/((x-.9)^2+.0