数值计算-第七讲课件

数值分析与算法(7)NumericalAnalysis&Algorithms清华大学计算机系喻文健数值分析与算法(7)NumericalAnalysisWenjianYu2第七章数值积分与数值微分

WenjianYu2第七章数值积分与数值微分

WenjianYu3数值积分的基本概念WenjianYu3数值积分的基本概念WenjianYu4数值积分目的与用途经典问题:算几何形体的面积、

体积，力学中

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

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

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

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

WenjianYu4数值积分目的与用途由一块平整的铝板压制WenjianYu5数值积分基本思想

...

积分系数

积分节点希望用较少的计算量得到较准确的结果WenjianYu5数值积分基本思想

...

积WenjianYu6插值型求积公式

中矩形公式梯形公式WenjianYu6插值型求积公式

中矩形WenjianYu7积分余项与代数精度

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

反映了计算的截断WenjianYu8积分余项与代数精度

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

求积公式WenjianYu8积分余项与代数精度

(至少0次代WenjianYu9求积公式的收敛性与稳定性

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

积分问题一般不太敏感

WenjianYu9求积公式的收敛性与稳定性

(一系列求积WenjianYu10求积公式的收敛性与稳定性

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

要尽量寻求稳定的求积公式WenjianYu10求积公式的收敛性与稳定性

这是控制WenjianYu11牛顿-柯特斯公式WenjianYu11牛顿-柯特斯公式WenjianYu12Newton-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系数表

一系列求积公式便于使用WenjianYu12Newton-Cotes公式

就是WenjianYu13Newton-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时的特例.例题(板书)WenjianYu13Newton-Cotes公式

n=1WenjianYu14Newton-Cotes公式

关键看积分:

(n阶公式至少有n次代数精度)一般不用n=3对应的N-C公式WenjianYu14Newton-Cotes公式

WenjianYu15低阶N-C公式的积分余项

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

2

详细过程自己看书WenjianYu15低阶N-C公式的积分余项

不保号WenjianYu16稳定性、收敛性

n=8,实际只使用n<8的偶数阶N-C公式,也看出代数精度不是越高越好WenjianYu16稳定性、收敛性

n=8,实际只WenjianYu17复合求积公式WenjianYu17复合求积公式WenjianYu18复合求积公式

(compositequadrature)

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

仍是机械求积公式WenjianYu18复合求积公式

(compositeWenjianYu19复合求积公式

2阶准确度

WenjianYu19复合求积公式

2阶准确度

WenjianYu20复合求积公式

与复合梯形公式对比,例7.4WenjianYu20复合求积公式

与复合梯形公式WenjianYu21复合求积公式步长折半的复合求积公式计算积分余项公式包含被积函数的高阶导数,很难应用.

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

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

只需再计算新增节点的函数值WenjianYu21复合求积公式步长折半的复合求积公式计WenjianYu22复合求积公式步长折半的复合求积公式计算复合Simpson公式的情况很少使用中矩形公式的原因

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

可类似构造复合中矩形公式,但在步长折半时,无法重用以前的结果WenjianYu22复合求积公式步长折半的复合求积公式计WenjianYu23Remberg积分算法WenjianYu23Remberg积分算法WenjianYu24复合梯形公式的余项展开式

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

所有小区间的积分求和:WenjianYu24复合梯形公式的余项展开式

(可与RWenjianYu25复合梯形公式的余项展开式

Th7.5

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

WenjianYu25复合梯形公式的余项展开式

Th7.5WenjianYu26RichardsonExtrapolation

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

WenjianYu26RichardsonExtrapoWenjianYu27RichardsonExtrapolation

WenjianYu27RichardsonExtrapoWenjianYu28Romberg算法龙贝格算法列三角形表格,按行依次计算计算公式

可证明:

具有2k+1次代数精度

(类似高阶差商的计算)WenjianYu28Romberg算法龙贝格算法

WenjianYu29Romberg算法

做等距的节点分布,充分利用其上函数值得到最高准确度的结果要求被积函数充分光滑!WenjianYu29Romberg算法

做等距的节点WenjianYu30Romberg算法

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.9460833WenjianYu30Romberg算法

011/80.9WenjianYu31Romberg算法

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

WenjianYu31Romberg算法

10.50000WenjianYu32自适应积分算法WenjianYu32自适应积分算法WenjianYu33自适应积分算法(adaptivequadrature)基本思想Romberg算法效果不好的情况积分节点没必要均匀分布怎样自动地非均匀取点,使计算

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

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

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

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

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

演示模块7.3,quadguiWenjianYu34自适应积分算法一个自适应求积算法

实WenjianYu35自适应积分算法一个自适应求积程序quadtx例子:更多讨论通过阈值设置控制相对误

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

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

[Q,fcnt]=quadtx(@humps,0,1,1e-3)fcnt=69,用了69个积分点WenjianYu35自适应积分算法一个自适应求积程序quWenjianYu36高斯求积公式WenjianYu36高斯求积公式WenjianYu37高斯求积公式

解得:

WenjianYu37高斯求积公式

解得:

WenjianYu38高斯求积公式

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

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

WenjianYu38高斯求积公式

高斯积分有2n+1次WenjianYu39高斯求积公式

Th7.7比自适应积分算法使用方便WenjianYu39高斯求积公式

Th7.7比WenjianYu40高斯-勒让德公式

012345高斯-勒让德积分表WenjianYu40高斯-勒让德公式

012345高斯-WenjianYu41高斯-勒让德公式

WenjianYu41高斯-勒让德公式

WenjianYu42数值微分WenjianYu42数值微分WenjianYu43数值微分

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

WenjianYu43数值微分

(向前差分)(WenjianYu44数值微分

WenjianYu44数值微分

WenjianYu45数值微分

h0.10.4516049081

0.050.45407616940.4548999231

00250.45469262880.45489811520.4548979947p359

准确的有效数字位数WenjianYu45数值微分

h0.10.45WenjianYu46数值微分的应用

二阶中心差分

例如n=2

WenjianYu46数值微分的应用

二阶中心WenjianYu47Matlab中的积分计算WenjianYu47Matlab中的积分计算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')函数句柄Matlab中的数值积分指定被积函数>>f1=inliMatlab中的数值积分一维积分的命令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-6Matlab中的数值积分一维积分的命令>>Q=quadMatlab中的数值积分二重、三重积分dblquad,quad2d,triplequad符号积分定义符号变量:

sym(),syms,符号积分:intsimple(表达式化简)double离散数据点积分复合梯形法trapz(x,y)>>symsx>>h=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)-6>>I=int(h)%不定积分I=10*atan(10*x-3)+5*atan(5*x-9/2)-6*x>>D=simple(int(h,0,1))%定积分D=5*atan(16/13)-6+10*pi>>Qexact=double(D)Qexact=29.858不用逐项运算符Matlab演示Matlab中的数值积分二重、三重积分>>symsx不用数值分析与算法(7)NumericalAnalysis&Algorithms清华大学计算机系喻文健数值分析与算法(7)NumericalAnalysisWenjianYu52第七章数值积分与数值微分

WenjianYu2第七章数值积分与数值微分

WenjianYu53数值积分的基本概念WenjianYu3数值积分的基本概念WenjianYu54数值积分目的与用途经典问题:算几何形体的面积、

体积，力学中

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

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

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

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

WenjianYu4数值积分目的与用途由一块平整的铝板压制WenjianYu55数值积分基本思想

...

积分系数

积分节点希望用较少的计算量得到较准确的结果WenjianYu5数值积分基本思想

...

积WenjianYu56插值型求积公式

中矩形公式梯形公式WenjianYu6插值型求积公式

中矩形WenjianYu57积分余项与代数精度

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

反映了计算的截断WenjianYu58积分余项与代数精度

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

求积公式WenjianYu8积分余项与代数精度

(至少0次代WenjianYu59求积公式的收敛性与稳定性

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

积分问题一般不太敏感

WenjianYu9求积公式的收敛性与稳定性

(一系列求积WenjianYu60求积公式的收敛性与稳定性

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

要尽量寻求稳定的求积公式WenjianYu10求积公式的收敛性与稳定性

这是控制WenjianYu61牛顿-柯特斯公式WenjianYu11牛顿-柯特斯公式WenjianYu62Newton-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系数表

一系列求积公式便于使用WenjianYu12Newton-Cotes公式

就是WenjianYu63Newton-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时的特例.例题(板书)WenjianYu13Newton-Cotes公式

n=1WenjianYu64Newton-Cotes公式

关键看积分:

(n阶公式至少有n次代数精度)一般不用n=3对应的N-C公式WenjianYu14Newton-Cotes公式

WenjianYu65低阶N-C公式的积分余项

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

2

详细过程自己看书WenjianYu15低阶N-C公式的积分余项

不保号WenjianYu66稳定性、收敛性

n=8,实际只使用n<8的偶数阶N-C公式,也看出代数精度不是越高越好WenjianYu16稳定性、收敛性

n=8,实际只WenjianYu67复合求积公式WenjianYu17复合求积公式WenjianYu68复合求积公式

(compositequadrature)

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

仍是机械求积公式WenjianYu18复合求积公式

(compositeWenjianYu69复合求积公式

2阶准确度

WenjianYu19复合求积公式

2阶准确度

WenjianYu70复合求积公式

与复合梯形公式对比,例7.4WenjianYu20复合求积公式

与复合梯形公式WenjianYu71复合求积公式步长折半的复合求积公式计算积分余项公式包含被积函数的高阶导数,很难应用.

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

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

只需再计算新增节点的函数值WenjianYu21复合求积公式步长折半的复合求积公式计WenjianYu72复合求积公式步长折半的复合求积公式计算复合Simpson公式的情况很少使用中矩形公式的原因

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

可类似构造复合中矩形公式,但在步长折半时,无法重用以前的结果WenjianYu22复合求积公式步长折半的复合求积公式计WenjianYu73Remberg积分算法WenjianYu23Remberg积分算法WenjianYu74复合梯形公式的余项展开式

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

所有小区间的积分求和:WenjianYu24复合梯形公式的余项展开式

(可与RWenjianYu75复合梯形公式的余项展开式

Th7.5

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

WenjianYu25复合梯形公式的余项展开式

Th7.5WenjianYu76RichardsonExtrapolation

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

WenjianYu26RichardsonExtrapoWenjianYu77RichardsonExtrapolation

WenjianYu27RichardsonExtrapoWenjianYu78Romberg算法龙贝格算法列三角形表格,按行依次计算计算公式

可证明:

具有2k+1次代数精度

(类似高阶差商的计算)WenjianYu28Romberg算法龙贝格算法

WenjianYu79Romberg算法

做等距的节点分布,充分利用其上函数值得到最高准确度的结果要求被积函数充分光滑!WenjianYu29Romberg算法

做等距的节点WenjianYu80Romberg算法

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.9460833WenjianYu30Romberg算法

011/80.9WenjianYu81Romberg算法

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

WenjianYu31Romberg算法

10.50000WenjianYu82自适应积分算法WenjianYu32自适应积分算法WenjianYu83自适应积分算法(adaptivequadrature)基本思想Romberg算法效果不好的情况积分节点没必要均匀分布怎样自动地非均匀取点,使计算

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

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

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

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

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

演示模块7.3,quadguiWenjianYu34自适应积分算法一个自适应求积算法

实WenjianYu85自适应积分算法一个自适应求积程序quadtx例子:更多讨论通过阈值设置控制相对误

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

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

[Q,fcnt]=quadtx(@humps,0,1,1e-3)fcnt=69,用了69个积分点WenjianYu35自适应积分算法一个自适应求积程序quWenjianYu86高斯求积公式WenjianYu36高斯求积公式WenjianYu87高斯求积公式

解得:

WenjianYu37高斯求积公式

解得:

WenjianYu88高斯求积公式

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

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

WenjianYu38高斯求积公式

高斯积分有2n+1次WenjianYu89高斯求积公式

Th7.7比自适应积分算法使用方便WenjianYu39高斯求积公式

Th7.7比WenjianYu90高斯-勒让德公式

012345高斯-勒让德积分表WenjianYu40高斯-勒让德公式

012345高斯-WenjianYu91高斯-勒让德公式

WenjianYu41高斯-勒让德公式

WenjianYu92数值微分WenjianYu42数值微分WenjianYu93数值微分

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

WenjianYu43数值微分

(向前差分)(WenjianYu94数值微分

WenjianYu44数值微分

WenjianYu95数值微分

h0.10.4516049081

0.050.45407616940.4548999231

00250.45469262880.45489811520.4548979947p359

准确的有效数字位数WenjianYu45数值微分

h0.10.45WenjianYu96数值微分的应用

二阶中心差分