实验6数值积分
数值分析第六章
可见， 积分系数出现了负数，而被积函数是恒正的函数，说明数值积分出现了不稳定现象.
6.1.3复化求积法
由定积分知识，定积分只与被积函数和积分区间 有关，而在对被积函数做插值逼近时，多项式的次 数越高，对被积函数的光滑程度要求也越高，且会 出现龙格(Runge)现象.如 时，牛顿-柯特斯公式 就是不稳定的. 因而，人们把目标转向积分区间， 类似分段插值，把积分区间分割成若干小区间，在 每个小区间上使用次数较低的牛顿-柯特斯公式， 然后把每个小区间上的结果加起来作为函数在整个 区间上积分的近似，这就是复化的基本思想.
j i
b a
f ( n1) ( ) b n f ( x)dx (b a) Ci f ( xi ) ( x xi )dx, (a, b). a ( n 1)! i 0 i 0
n
(6-10) 称公式(6-10)为牛顿-柯特斯公式的截断误差 . 柯特斯系数 与被积函数和积分区间都无关， 只要给出区间等分数 即可求出.
Matlab程序与计算实例
计算步骤： (a) 计算柯特斯系数中的积分表达式
li
j i
t j , (i 0, 1, 2, , n). j 0 i j
n
(b) 计算柯特斯系数
j i
n n (1)ni Ci (t j )dt. 0 n i !(n i )! j 0
6.1 实验的主要方法
插值型求积法 复化求积法 龙贝格（Romberg）求积法 高斯(Gauss)求积法 二重积分的常用数值积分方法简介
6.1.1数值积分的基本思想
数值积分法的思想来源于定积分的定义，考虑带权 定积分 ( x) f ( x)dx lim ( x ) f ( x )x . ( x) 0 其中 max{x }， 为权函数.当 ( x) 1 时，即是普通积 分. 用f ( x) 在点 a x x x b 处的函数值 f ( x ) (i 0, 1, , n) 的线 性组合 I ( f ) A f ( x ) A f ( x ) A f ( x ) A f ( x ). 作为积分I ( f ) ( x) f ( x)dx 的近似值，即 I ( f ) ( x ) f ( x )dx A f ( x ). (6-1) 称(6-1)为求积公式，称 R( f ) I ( f ) I ( f ). (6-2) x A (i 0, 1, 2, , n) 分别 为求积公式(6-1)的余项或误差，及 称为求积公式(6-1)的求积节点和求积系数. 这里 求积系数 A (i 0, 1, 2, , n) 只与权函数 ( x) 及积分区间[a, b] 有 关，而与 f ( x) 无关.这类求积方法通常称为机械求积 法。
6.数值积分与数值微分
b Ak =∫a lk( x)dx
2、n=2时 有(x0，f（x0）)， (x1，f（x1）) ， (x2，f（x2）) b ∫a f( x)dx ≈ A0 f（x0）+ A1f（x1）+ A2f（x2） = A0 f（a）+ A1f（x1）+ A2f（b） b b A0 =∫a l0 ( x)dx =∫a [(x- x1)( x- x2)/[(x0-x1)(x0-x2)]] dx b =∫a [(x-x1)( x-b)/[(a-x1)(a-b)]]dx =[(b-a)(2a+b-3x1)]/[6(a-x1)]
b b
3. n=4时,有5个节点 4 （ 4 ） 4 C0 =(-1) /[4(0!)(4!)] ∫0 (t-1)(t-2)(t-3)(t-4)dt=7/90 4 C1（4） =(-1) 3 /[4(1!)(3!)] ∫0 t(t-2)(t-3)(t-4)dt=32/90 4 C2（4） =(-1) 2 /[4(2!)(2!)] ∫0 t(t-1)(t-3)(t-4)dt=12/90 4 （ 4 ） 1 C3 =(-1) /[4(3!)(1!)] ∫0 t(t-1)(t-2)(t-4)dt=32/90 4 C4（4） =(-1) 0 /[4(4!)(0!)] ∫0 t(t-1)(t-2)(t-3)dt=7/90 可以看出 C0（4） =C4(4) ， C1(4) = C3(4) ∴ I4 =(b-a)∑ C k (4) f(x0 +kh) =[(b-a)/90][7f(x0)+32 f(x1)+12f(x2)+32f(x3)+7(x4)] 这是Cotes求积公式。 Cotes系数性质： b (1)与[a,b]无关，f(x)=1时， (b-a)∑Ck（n）1=∫a 1 dx=b-a 所以 ∑Ck（n） = 1 (2)对称性Ck（n） =Cn-k(n)
数值计算方法之数值积分
数值计算方法之数值积分数值积分是数值计算中的一个重要内容，它是对函数在其中一区间上的积分进行数值近似计算的方法。
数值积分在计算机科学、自然科学以及工程领域都有广泛的应用，如求解不定积分、概率密度函数的积分、求解微分方程初值问题等。
数值积分的基本思想是将积分区间划分为若干小区间，然后对每个小区间进行数值近似计算，再将结果相加得到近似的积分值。
常用的数值积分方法包括矩形法、梯形法、辛普森法等。
首先介绍矩形法。
矩形法是将积分区间划分为若干个小区间，然后用每个小区间的函数值与该小区间的宽度相乘得到每个小矩形的面积，最后将所有小矩形的面积相加得到近似的积分值。
矩形法分为左矩形法、右矩形法和中矩形法三种。
左矩形法即用每个小区间的最左端点的函数值进行计算，右矩形法用最右端点的函数值进行计算，中矩形法用每个小区间中点的函数值进行计算。
梯形法是将积分区间划分为若干个小区间，然后用每个小区间两个端点的函数值与该小区间的宽度相乘，再将每个小梯形的面积相加得到近似的积分值。
梯形法相较于矩形法更为精确，但需要更多的计算量。
辛普森法是将积分区间划分为若干个小区间，然后用每个小区间的三个点的函数值进行插值，将插值函数进行积分得到该小区间的近似积分值，最后将所有小区间的近似积分值相加得到近似的积分值。
辛普森法相比矩形法和梯形法更为精确，但计算量更大。
除了以上几种基本的数值积分方法外，还有龙贝格积分法、高斯积分法等更为精确的数值积分方法。
这些方法的原理和步骤略有不同，但都是通过将积分区间分割为若干小区间，然后进行数值近似计算得到积分值的。
总结起来，数值积分是通过将积分区间分割为若干小区间，然后对每个小区间进行数值近似计算得到积分值的方法。
常用的数值积分方法包括矩形法、梯形法、辛普森法等。
数值积分在计算机科学、自然科学以及工程领域均有广泛应用，是数值计算中的重要内容。
数值分析第6章积分
其中(-1,1).
b
ò f (x )dx
a
»
å j
n
=0
A j f (x j )
(6.1)
定理6.1 求积公式(6.1)为插值型求积公式的 充要条件是它的代数精度至少为n次. 证:先证必要性 设(6.1)是插值型的,则
b
R n [f ] =
ò
a
f
(x ) wn + 1 (x )dx (n + 1)!
顿-柯特斯（Newton-Cotes）求积公式. 下面推导N-C求积公式的求积系数公式.
根据求积系数计算公式(6.4)有
1 Aj = 蝌 l j ( x)dx = w¢ (x j ) a
b b a
wn+ 1 ( x) dx x- xj
令积分变换 x=a + t h, 则
wn+ 1 ( x) = h n+ 1t (t - 1) L (t - n),
òl
a
j
( x) dx
(j=0,1,2, L ,n)
(6.4)
若求积公式(6.1)中的求积系数具有(6.4)的形 式,则称(6.1)为插值型求积公式.
插值型求积公式(6.3)的截断误差为
b b
Rn [ f ] =
R ( x)dx = 蝌
n a a
f ( n+ 1) (x ) wn+ 1 ( x)dx (6.5) (n + 1)!
于是
Aj =
n- j n ¢ wn ( x ) = ( 1) j !( n j )! h , +1 j
(- 1) h t (t - 1) L (t - j + 1)(t - j - 1) L (t - n)dt ò j !(n - j )! 0
第六章 数值积分
由(5)式，显而易见+1) (ε ) = 0 , 可知，R(1, f)=0，所以我们所n+1点的求积公式(3)至少具 有n次的代数精确度．进一步可以证明，当n为偶数时， 求积公式(3)的代数精确度可以达到n+1次．
三、几种常见的Newton-Cotes求积公式 对 n=0, 1, 2, 按公式(3)可以得出下面三种常见的Newton-Cotes 求积公式. 1. n=0时的矩形求积公式 分别以积分区间[a, b]的左、右端点和区间中点，即x=a，b， (a+b)/2为求积节点得到： 左矩形求积公式： I ( f ) = f (a)(b − a) + R 1 右矩形求积公式： 中矩形求积公式： I ( f ) = f ( a + b)(b − a) + R 3 2 三个求积公式的误差估计，可将函数f(x) 分别在 x = a, b, a +b 处展开到含f(x)的一阶导数的Taylor公式在区间[a, b]上积分 推得．
b−a a +b { f (a) + 4 f ( ) + f (b)}+ R(σ, f ) 6 2 b5 − a5 （1）当 f (x) = x4时，I ( f ) = 54 b − a 4 (a +b) I3( f ) = (a + +b4 ) ≠ I ( f ) 6 4 I ( f ) = I3( f ) + R(σ, f ) =
定义１ 若对任意的 ， pn (x) ∈P [a, b] n 求积公式(2)的误差都满足 ，则称 该求积公式具有n次代数精确度． 验证一个求积公式所具有的代数精确度用 定义１是极不方便的，为此给出另一个定义．
定义2 而
若对函数 ，
数学实验之数值积分PPT文档共19页
数学实验之数值积分
11、用道德的示范来造就一个人，显然比用法律来约束他更有价值。—— 希腊
12、法律是无私的，对谁都一视同仁。在每件事上，她都不徇私情。—— 托马斯
13、公正的法律限制不了好的自由，因为好人不会去做法律不允许的事 情。——弗劳德
14、法律是为了保护无辜而制定的。——爱略特 15、像房子一样，法律和掉进坑里。——黑格尔 32、希望的灯一旦熄灭，生活刹那间变成了一片黑暗。——普列姆昌德 33、希望是人生的乳母。——科策布 34、形成天才的决定因素应该是勤奋。——郭沫若 35、学到很多东西的诀窍，就是一下子不要学很多。——洛克
数值计算方法 第4版 第6章 数值积分和数值微分
m 次的代数精度。
定义 求积公式
对 f (x) 1, x, x2,
ab
f
( x)dx
n
Ak
f
(xk )
k 0
, xm ，准确成立，而对于 f (x) xm1 ，不能准确成立，
则称该求积公式
ab
f
( x)dx
n
Ak
f
(xk )
k 0
具有 m 次的代数精度。
求积公式
ab
f
( x)dx
n
Ak
, xn
，使求积公式 ab
f
( x)dx
n
Ak
f
(xk )
准
k 0
A0 A1 An b a
A0
x0
A1x1
An xn
1 2
(b 2
a2)
A0
x0n
A1 x1n
An xnn
1 (bn1 n 1
an1)
这是关于Ak 的线性方程组，其系数矩阵是凡德蒙矩阵，当节点xk 互
异时， Ak , k 0,1, , n 有唯一解。即对 f (x) 1, x, x2 , , xn 求积公式
定理 对于给定的 n+1 个节点 xk , k 0,1, , n 总存在求积系数
Ak , k 0,1,
精度。
, n 使求积公式 ab
n
f (x)dx
Ak
f (xk ) 至少有
n
次的代数
k 0
证明此时 Ak , k 0,1, , n 有唯一解即可。
证 令 f (x) 1, x, x2 ,
确成立，即
a
Ak f (xk )
k 0
6 数值积分与微分
n
( ) f x （ f xi）lin1 ( x) n ( x) x [a, b] i 1 n! (n) n b b b f ( ) f x dx （ f xi） lin1 ( x)dx n ( x)dx a a a i 1 n!
a
ba f x dx f a f b 2
常用的Newton-Cotes公式 梯形公式
b
a
ba f x dx f a f b 2
Simpson公式或抛物线公式
b
a
ba f x dx f a 4 f 6
2
1
可以事先算出cotes系数，如n=2的两个系数为
1 t 2 1 1 C1 dt 1 t dt 0 1 2 0 2
2
1
C2
1 t 1 1 1 dt t dt 0 2 1 0 2
所以2 点的Newton-Cotes公式为
©
b
的代数精度。
取f x x R( x )
2 2
2
1
1 1 2 1 2 x dx 1 2 0 2 6 2
2
故本题求积公式代数精度为1。
©
例 6.2 确定求积公式
2h
2 h
f ( x)dx Af (h) Bf (0) Cf (2h)
的参数A，B，C，使它具有尽可能高的代数精 度，并指出相应的代数精度。
解：要先求出具体的求积公式，再判断代数精度。
依次取f ( x) 1, x, x2代入求积公式的两端，并令其相等
A B C 4h A 2C 0 A 4C 16h / 3
数值计算06-数值积分与数值微分
用 inline 函数定义被积函数： >> f=inline('2/sqrt(pi)*exp(-x.^2)','x'); >> y=quad(f,0,1.5)
y= 0.9661
• 矩形区域上的二重积分的数值计算
I yM xM f (x, y)dxdy ym xm
格式： 矩形区域的双重积分： y=dblquad(Fun,xm,xM,ym,yM)
数值计算
第六章 数值积分与数值微分
1
§6.1 引 言
一、数值积分的必要性
讨论如下形式的一元函数积分
b
I ( f ) f (x)dx
a
在微积分里，按Newton-Leibniz公式求定积分
b
I ( f ) a f (x)dx F (b) F (a)
要求被积函数 F x
☞ 有解析表达式；
☞ f x的原函数 F x 为初等函数．
k 0
称为求积公式 余项(误差).
构造或确定一个求积公式，要解决的问题包括:
(i) 确定求积系数 Ak 和求积节点 xk；
(ii) 确定衡量求积公式好坏的标准； (iii) 求积公式的误差估计和收敛性分析.
数值积分的基本问题
针对某类函数，选择合适的求积结点和求积系数，使得求积 公式(1) 具有尽可能小的截断误差或尽可能高的代数精度。
2
若f ( x)在区间[a,b]上有四阶连续导数。则Simpson
公式的截断误差
R2
(b a)5
2880
f (4)( ) (a,b)
(6.3.8)
且具有三次代数精度。
Simpson3/8公式：
数 值 积 分 法
Fct(t,h,x,f); k1[0]=f[0]; k1[1]=f[1]; x1[0]=x[0]+h/2*k1[0]; x1[1]=x[1]+h/2*k1[1]; Fct(t+h/2,h,x1,f); k2[0]=f[0]; k2[1]=f[1]; x1[0]=x[0]+h/2*k2[0]; x1[1]=x[1]+h/2*k2[1]; Fct(t+h/2,h,x1,f); k3[0]=f[0]; k3[1]=f[1]; x1[0]=x[0]+h*k3[0]; x1[1]=x[1]+h*k3[1]; Fct(t+h,h,x1,f); k4[0]=f[0]; k4[1]=f[1]; x1[0]=x[0]+h/6*(k1[0]+2*k2[0]+2*k3[0]+k4[0]); x1[1]=x[1]+h/6*(k1[1]+2*k2[1]+2*k3[1]+k4[1]); x[0]=x1[0]; x[1]=x1[1]; } void Fct(double t, double h, double* x, double* f) { double g = 9.81; double den; f[0] = x[1]; den = lagrange(x); f[1] = -5*g+0.5*0.02*3.14*0.05*0.05*x[1]*x[1]*den; } double lagrange(double* x) {
五 实验结果及分析
1.实验结果图表 欧拉法
二阶龙哥库塔法
四阶龙哥库塔法
2.在程序循环到高度小于零时速度值已不可信， 这是因为在编程时未考虑到高度 小于零时的公式符号变化，此时数据已无实际意义，所以我没有修改； 3.小球落地时的速度值可由 txt 文件查询到分别为 欧拉法 354.883886m/s 二阶龙哥库塔法 354.869480m/s 四阶龙哥库塔法 354.844668m/s 小球落地时循环步数可由 txt 文件查询到，再乘以步长即可得到时间 欧拉法 59.22s 二阶龙哥库塔法 59.22s 四阶龙哥库塔法 59.21s 均为落地前一秒的数值近似得到
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
1
实验六 数值积分(Quadrature)
6.1 实验目的
掌握数值积分的代数精度与收敛阶的含义，会利用matlab求解符号积分和
数值积分，并会借助数学软件Matlab求解一些简单的实际问题。
6.2 实验内容
1、Matlab中求解符号积分和数值积分的方法；
2、建立飞船的轨道周长和手的面积等实际问题的数学模型，并借助数学软
件Matlab求解．
6.3 实验步骤
6.3.1 函数表达式已知时的积分
Matlab中求积分的方法包括符号积分（精确计算）和数值积分两种。前者
只有一个命令int（注：求导数的命令为diff）,后者命令较多，详见表1。
表1 Matlab中数值积分函数一览表
函数名 调用格式 方法说明
适应范围（详
见Help）
quad q=quad(fun,a,b,tol) 自适应Simpson 法 被积函数不光滑，低精度要
求
quadl q=quadl(fun,a,b,tol) 自适应 Lobatto 法 被积函数光滑，高精度要
求
quadgk
[q,errbnd]=
quadgk(fun,a,b,tol)
自适应Gauss-
Kronrod 法
震荡型被积函
数，高精度要
求，积分限可
为无穷
trapz q=trapz(X,Y) 梯形法
离散数据点积
分
dblquad
q = dblquad(fun,
xmin,xmax,ymin,ymax)
调用一维数值积分方法 矩形区域上的
二重积分
quad2d q = quad2d(fun,a,b,c,d)
平面区域上的
二重积分
triplequad
triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax) 调用一维数值积分方法 三维长方形区
域上的三重积
分
2
教师示范1：符号积分
int是matlab中符号工具箱的一个命令，它能精确计算不定积分和定积分，
微积分教材中的积分一般均可用该种方法积出。
例1 求不定积分222(1)xdxx。
Matlab命令序列为：
syms x
s=int(-2*x/(1+x^2)^2)
运行结果为：1/(x^2 + 1)，即
222
21(1)1xdxxx
。
验算所得结果，即验证
222
121(1)dxdxxx
，
Matlab命令为：diff(s)，运行结果为：-(2*x)/(x^2 + 1)^2.
其余见文件：integrate_symbolic.m。
教师示范2：数值积分
例2 求不定积分21.9620xedx
Matlab命令为：
format long
q1 = quad(@(x) exp(-x.^2/2),0,2)
q2 = quadl(@(x) exp(-x.^2/2),0,2)
q3 = quadgk(@(x) exp(-x.^2/2),0,2)
结果分别为1.196288044999181，1.196288103140270和1.196288013322608。
实际上，
2
1.96
2
0
1.19628801332260820293142377047xedx
由此，函数quadgk比quad、quadl的计算精度要高。
6.3.2 离散数据的数值积分
当函数的表达式未知，但已知函数在一系列节点处的函数值时，()bafxdx的
计算有三种方法。方法1，利用表1中的梯形法；方法2，先用样条插值函数逼
3
近未知函数，然后利用表1中的命令求该样条函数的积分；方法3，如果所求的
图形可视为以点(,)iixy为顶点的多边形，则可利用多边形的面积公式（1）
11111112nkkkknnkAxyxyxyxy
（1）
直接求解。
具体见例3。
例3 用不同方法计算不定积分100cossin10xdx，比较不同方法的精度。
程序如下。
%% 比较不同方法求定积分的精度
a = 0; b = 10;
syms x;
exact = vpa( int(cos(x),a,b), 20); % 符号定积分，小数点20位的精确值
q = zeros(7,1);
q(1) = quad(@cos,a,b);
q(2) = quadl(@cos,a,b);
q(3) = quadgk(@cos,a,b);
% 离散数据后，用梯形法求积分
t = a:0.1:b;
y = cos(t);
q(4) = trapz(t, y);
% 用三次样条近似函数后，求积分
pp = spline( t, y );
q(5) = quad(@(t)ppval(pp,t),a,b);
q(6) = quadl(@(t)ppval(pp,t),a,b);
q(7) = quadgk(@(t)ppval(pp,t),a,b);
%
disp('精确值')
disp(exact)
format long;
disp('近似值')
disp(q);
disp('各种方法的误差')
for i = 1:7
fprintf(1,'%12.10e\n',double( q(i)-exact ) );
end
以上程序运行结果见表2。
4
表2 Matlab中不同数值积分方法精度比较
函数 结果 方法说明 误差
int -0.5440211108893698134
符号积分，精
确解
quad -0.544021121312844 自适应Simpson 法 -1.04e-008
quadl -0.544021110988547 自适应 Lobatto 法 -9.92e-011
quadgk -0.544021110889370 自适应Gauss- Kronrod 法 -1.83e-016
trapz -0.543567684387146 梯形法 4.53e-004
spline quad -0.544021095007636 三次样条近似 1.59e-008 quadl -0.544020920620326 1.91e-007
quadgk -0.544020997819316 1.13e-007
详见文件：integrate_numeric.m。
6.3.2 学生练习
学生练习1：网球拍的拍面外形可近似看成椭圆，如图1所示。试计算该网
球拍拍面外形的周长。首先给出求解椭圆周长的积分公式，然后数值求解。
图1. 网球拍的拍面
据中国载人航天工程网报道[1]，中国第一艘载人航天飞船神舟五号的初始轨
道为近地点200公里、远地点350公里的椭圆轨道，试计算该椭圆轨道的长度。
设地球为半径6371公里的球形。
（选做）你可以通过查阅权威网站(如中国航天网)提供的数据，计算嫦娥
二号绕月飞行的椭圆轨道的长度。
5
学生练习2：实验5插值法的学生实验练习2中已绘制出你的手，试计算你
的手的面积。
图2. 一只手
根据文献[2]，“用病人的手估计估计烧伤面积已成为诊断烧伤面积的公认方
法之一”，并指出“男性手掌的面积大致占人体表面积的0.81%，女性则大致为
0.67%”。试由此估计你的体表面积。
提示：手的面积的求法有多种。一种最简单的方法是利用多边形的面积。可
选的方法是定积分或二重积分，请同学们认真思考。给出不少于两种的计算方法。
6.4 思考题（选做）
1．用数值积分方法计算含参变量的积分
1110(,)(1)zwzwttdt，10()xtxtedt
，
并将你的计算结果与Matlab中的内置函数beta(z,w)和gamma(x)的结果进行比
较。前者可参见integrate_numeric.m。
2．查阅资料学习如何用Monte-Carlo算法计算二重积分。然后，分别用不
同的100,500,1000,...N，计算积分
xyDIedxdy
，其中[1,2][2,3]D，
对于每个N计算l次，计算出l次结果的均方差()N，做出()NN之间的曲
线。
参考文献
[1] 神舟五号载人飞船简介[DB]. http://www.cmse.gov.cn/project/show.php?itemid=
114，2008-09-17.
[2] 张向清. 一只手的面积有多大[J]. 中国烧伤创疡杂志，1997,2：23.