1 、牛顿 - 莱布尼兹公式 另外若给出的函数 f(x) 是数据表，也不好求函数的积分。 计算定积分的方法： 但是求函数 f(x) 的原函数 F(x) 不一定比计算积分容易， 例如函数 找不到用初等函数表示的原函数。 一、数值求积的基本思想 实验 4 数值积分与微分 主讲人：魏志强.

Presentation on theme: "1 、牛顿 - 莱布尼兹公式 另外若给出的函数 f(x) 是数据表，也不好求函数的积分。 计算定积分的方法： 但是求函数 f(x) 的原函数 F(x) 不一定比计算积分容易， 例如函数 找不到用初等函数表示的原函数。 一、数值求积的基本思想 实验 4 数值积分与微分 主讲人：魏志强."— Presentation transcript:

1

2 1 、牛顿 - 莱布尼兹公式 另外若给出的函数 f(x) 是数据表，也不好求函数的积分。 计算定积分的方法： 但是求函数 f(x) 的原函数 F(x) 不一定比计算积分容易， 例如函数 找不到用初等函数表示的原函数。 一、数值求积的基本思想 实验 4 数值积分与微分 主讲人：魏志强

3 数学实验课件

4 2 、积分中值定理 但是点 的具体位置一般不知道，故难以准确算出 的值。 有两种近似方法： 得到 梯形公式 另一种是矩形法， 左矩形公式 一种是梯形法，用 代替 中矩形公式 右矩形公式 简称矩形公式

5 二、 数值积分基本概念 定义 1 ：公式 叫做数 值求积公式。 其中 x k 称为求积节点， A k 称为求积系数, 亦称伴随节点 x k 的权。 定义 2 ：代数精度 若求积公式（ 1 ）对于次数不超过 m 的多项式均能准确成立，但对于 m+1 次的多项式 就不能准确成立，则称此求积公式的代数精度为 m 。 如何判断代数精度？

6 三、插值型的求积公式 设给定一组节点 : 且已知函数在这些节点处的函数值 f(x i )(i=0,1,…n), 由第二 章知可以作插值函数 L n (x), 由于 L n (x) 为多项式，所以其积 分可以很容易求得： 记 则上式 = 定义：称公式 其中 为求积系数。 为插值型的求积公式。 其余项为：

7 对等分的情况 ( k=0,1…,n ) 代入插值求积有 称为牛顿 - 柯特斯求积公式,C k （ n) 称为柯特斯系数 引进记号 ( k=0,1…,n ) 则

8 n = 1: 梯形公式 代数精度 = 1 n = 2: 辛普森公式 代数精度 = 3 余项

9 n = 4: 柯特斯公式 代数精度 = 5 余项

10 变步长求积方法 复化辛普森公式及其误差 将积分区间 [a,b] 划分为 n 等分, 记子区间 的中点为 在每个小区间上应用辛卜生 公式，则有 记 称为复化辛普森公式 ，余项

11 龙贝格求积公式 梯形法的递推化 为了提高求积精度，可在复化求积的基础上将积分小区间 [x k, x k+1 ] 二分一次，增加了一个分点 x k+1/2 = (x k +x k+1 ) /2, 记 h=(b-a)/n ，则： 称之为梯形法的递推化

12 龙贝格算法 由复化梯形公式的余项知： 有 即： 所以说明复化梯形公式二分前后两次计算值的线性组合 就为复化辛普森求积公式。

13 同理对辛普生公式进行二分处理，前后两次计算值的误差 进行比较， 复化柯特斯公式 有 即：

14 重复同样的操作，可进一步得到龙贝格（ Romberg) 公式： 高斯求积公式？

15 四、 数值积分的实现方法 1 ．变步长辛普森法 基于变步长辛普森法， MATLAB 给出了 quad 函数来求定积 分。该函数的调用格式为： [I,n]=quad('fname',a,b,tol,trace) 其中 fname 是被积函数名。 a 和 b 分别是定积分的下限和上限。 tol 用来控制积分精度， 缺省时取 tol=0.001 。 trace 控制是否展现积分过程，若取 非 0 则展现积分过程，取 0 则不展现，缺省时取 trace=0 。 返回参数 I 即定积分值， n 为被积函数的调用次数。

16 例 1 求定积分。 (1) 建立被积函数文件 fesin.m 。 function f=fesin(x) f=exp(-0.5*x).*sin(x+pi/6); (2) 调用数值积分函数 quad 求定积分。 [S,n]=quad('fesin',0,3*pi) S = 0.9008 n = 77 To Matlab

17 2 ．牛顿－柯特斯法 MATLAB 给出 quadl 函数来求定积分。 该函数的调用格式为： [I,n]=quadl('fname',a,b,tol,trace) 其中参数的含义和 quad 函数相似，只是 tol 的缺省值 取 10^(-6) 。该函数可以更精确地求出定积分的 值，且一般情况下函数调用的步数明显小于 quad 函数，从而保证能以更高的效率求出所需的定积 分值。 如： [S,n]=quadl('fesin',0,3*pi)

18 例 2 求定积分。 (1) 被积函数文件 fx1.m 。 function f=fx1(x) f=x.*sin(x)./(1+cos(x).*cos(x)); (2) 调用函数 quadl 求定积分。 I=quadl('fx1',0,pi) I = 2.4674

19 调用函数 quad 求定积分 ( 可以文件存储）： format long; fx=inline('exp(-x)'); [I,n]=quad(fx,1,2.5,1e-10) I = 0.28579444254766 n = 65 例 3 分别用 quad 函数和 quad8 函数求定积分的近 似值，并在相同的积分精度下，比较函数的调用 次数。

20 调用函数 quadl 求定积分： format long; fx=inline('exp(-x)'); [I,n]=quadl(fx,1,2.5,1e-10) I = 0.28579444254754 n = 33

21 4 ．被积函数是由表格定义时 在 MATLAB 中，对由表格形式定义的函数关系的求定积 分问题用 trapz(X,Y) 函数。其中向量 X,Y 定义函数关系 Y=f(X) 。 例 4 用 trapz 函数计算定积分。 命令如下： X=1:0.01:2.5; Y=exp(-X); % 生成函数关系数据向量 trapz(X,Y) ans = 0.28579682416393

22 五、二重定积分的数值求解 使用 MATLAB 提供的 dblquad 函数就可以直接求 出二重定积分的数值解。该函数的调用格式为： I=dblquad(f,a,b,c,d,tol,trace) 该函数求 f(x,y) 在 [a,b]×[c,d] 区域上的二重定积分。 参数 tol ， trace 的用法与函数 quad 完全相同。

23 例 5 计算二重定积分 (1) 建立一个函数文件 fxy.m ： function f=fxy(x,y) global ki; ki=ki+1; %ki 用于统计被积函数的调用次数 f=exp(-x.^2/2).*sin(x.^2+y); (2) 调用 dblquad 函数求解。 global ki;ki=0; I=dblquad('fxy',-2,2,-1,1) ki I = 1.57449318974494 ki = 1038

24 数值微分 一、 差分与差商 高阶差分、差商的定义

25 二、插值型求导公式 用插值多项式来构造数值微分的基本方法： 对给定的 f(x) 的函数表，构造对应的插值多项式 P n (x), 再令 去求函数微商。 常用的是在节点处带余项的数值微分公式： 由于 分析 : 故

26 1 、两点公式 设两个节点 x 0, x 1 处的函数值为： f(x 0 ), f(x 1 ), 求 f ′ (x 0 ), f ′ (x 1 ) 。 线性插值公式为 记 h= x 1 -x 0 有 余项：

27 2 、三点公式 设已知三点 x 0,x 1 =x 0 +h,x 2 =x 0 +2h 上的函数值 f(x 0 ), f(x 1 ), f(x 2 ), 求 f ′ (x 0 ), f ′ (x 1 ), f ′ (x 2 ) 。 对 t 求导可得： 令 由 故 令 t=0,1,2 ，便得到三节点处的导数： 余项：

28 二、 数值微分的实现 在 MATLAB 中，没有直接提供求数值导数的函数， 只有计算向前差分的函数 diff ，其调用格式为： DX=diff(X) ： 计算向量 X 的向前差分， DX(i)=X(i+1)-X(i) ， i=1,2,…,n-1 。 DX=diff(X,n) ： 计算 X 的 n 阶向前差分。例如， diff(X,2)=diff(diff(X)) 。 DX=diff(A,n,dim) ：计算矩阵 A 的 n 阶差分， dim=1 时 ( 缺省状态 ) ，按列计算差分； dim=2 ，按行 计算差分。

29 例 6 生成以向量 V=[1,2,3,4,5,6] 为基础的范得蒙矩阵，按列进 行差分运算。 命令如下： V=vander(1:6) DV=diff(V) % 计算 V 的一阶差分

30 例 7 用不同的方法求函数 f(x) 的数值导数，并在同一个坐标 系中做出 f'(x) 的图像。 程序如下： f=inline('sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2'); g=inline('(3*x.^2+4*x-1)./sqrt(x.^3+2*x.^2- x+12)/2+1/6./(x+5).^(5/6)+5'); x=-3:0.01:3; p=polyfit(x,f(x),5); % 用 5 次多项式 p 拟合 f(x) dp=polyder(p); % 对拟合多项式 p 求导数 dp dpx=polyval(dp,x); % 求 dp 在假设点的函数值 dx=diff(f([x,3.01]))/0.01; % 直接对 f(x) 求数值导数 gx=g(x); % 求函数 f 的导函数 g 在假设点的导数 plot(x,dpx,'-.',x,dx,'.',x,gx, '-' ); % 作图