第7章MATLAB数值积分与微分
7.1 数值积分
7.2 数值微分
7.1 数值积分
7.1.1 数值积分基本原理
求解定积分的数值方法多种多样，如简单的梯形法、辛普生 (Simpson)法、牛顿－柯特斯(Newton-Cotes)法等都是经常采用的方法。
它们的基本思想都是将整个积分区间[a,b]分成n 个子区间[x i ,x i+1]，i=1,2,…,n ，其中x 1=a ，x n+1=b 。这样求定积分问题就分解为求和问题。
7.1.2 数值积分的实现方法(Quadrature) 1．自适应辛普生法
基于变步长辛普生法，MATLAB给出了quad函数来求定积分。该函数的调用格式为：
[I,n]=quad(fun,a,b,tol,trace)
fun是被积函数名、句柄或内联函数对象;
a和b分别是定积分的下限和上限;
tol用来控制积分精度，缺省时取tol=10-6;
trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，缺省时取trace=0
I即定积分值，n为被积函数的调用次数。
例7-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
2．牛顿－柯特斯法
基于牛顿－柯特斯法，MATLAB给出了quad8函数来求定积分。该函数的调用格式为：[I,n]=quad8(fun,a,b,tol,trace)
其中参数的含义和quad函数相似。 该函数可以更精确地求出定积分的值，且一般情况下函数调用的步数明显小于quad函数，从而保证能以更高的效率求出所需的定积分值。
例7-2 求定积分。
(1) 被积函数文件fx.m。
f=inline('x.*sin(x)./(1+cos(x).*cos(x)) '); I=quad8(f,0,pi)
I =
2.4674
2．高阶方法
MATLAB给出了quadl函数来替代quad8求定积分[I,n]=quadl(fun,a,b,tol,trace)
其中参数的含义和quad函数相似。 该函数可以更精确地求出定积分的值，且一般情况下函数调用的步数明显小于quad函数，从而保证能以更高的效率求出所需的定积分值。
例7-3 分别用quad函数、quad8和quadl函数求定积分的近似值，并在相同的积分精度下，比较函数的调用次数。
调用函数quad求定积分：
format long;
fx=inline('exp(-x)');
[I,n]=quad(fx,1,2.5,1e-10)
I =
0.28579444254766 n =
65调用函数quadl求定积分：[I,n]=quadl(fx,1,2.5,1e-10)
I =
0.28579444254754
n =
33
调用函数quad8求定积分：[I,n]=quad8(fx,1,2.5,1e-10)
I =
0.28579444254754
n =
33
7.1.3 二重定积分的数值求解
二重积分可由单重积分函数来构建。
1、自己利用单重积分函数编程进行二重积分计算；
2、使用MATLAB提供的dblquad函数求二重定积分的数值解。
I=dblquad(f,a,b,c,d,tol,method)
该函数求f(x,y)在[a,b]×[c,d]区域上的二重定积分。参数tol的用法与函数quad完全相同。method 用于指定用哪个数值积分函数，如
@quad8,@quadl或用户自己定义的积分函数
例7-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
？？
已知二自变量函数在矩形网格下采样点的值，如何计算其积分？
x=-1:0.1:1;
y=-2:0.1:2;
[X,Y]=meshgrid(x,y);
Z=sqrt(X.^2+Y.^2);
mesh(X,Y,Z)
I1=trapz(y,trapz(x,Z,2))
I2=trapz(x,trapz(y,Z),2)
7.1.4 三重定积分的数值求解
三重积分可由单重积分函数来构建。
1、自己利用单重和二重积分函数编程进行三重积分计算；
2、使用MATLAB提供的triplequad函数求三重定积分的数值解。
triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol,method)
e.g.
Q = triplequad(inline('y*sin(x)+z*cos(x)'),0,pi,0,1,-1,1)
或
Q = triplequad(@integrnd,0,pi,0,1,-1,1)
where integrnd.m is the M-file
function f = integrnd(x,y,z)
f = y*sin(x)+z*cos(x);
一重积分计算平面图形面积，二重积分计算体积，三重积分的物理意义是什么？？如何使用triplequad计算一个非立方体的质量?
？？
已知三自变量函数在空间立方网格下采样点的值，如何计算其积分？x=-1:0.1:1;
y=-2:0.1:2;
z=1:0.1:2;
[X,Y,Z]=meshgrid(x,y,z);
F=sqrt(X.^2+Y.^2+Z.^2);
slice(X,Y,Z,F, [-0.5 0 0.5],[-1 0 1], 1.5)
shading interp
daspect([1 1 1])
camlight;
I1=trapz(z,trapz(y,trapz(x,Z,2)),3)
I2=trapz(z,trapz(x,trapz(y,Z),2),3)
I3=trapz(x,trapz(y,trapz(z,Z,3)),2)
例6-6 生成以向量V=[1,2,3,4,5,6]为基础的范得蒙矩阵，按列进行差分运算。
命令如下：
V=vander(1:6)
DV=diff(V) %计算V的一阶差分
例6-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, 'r-'); %作图
标量场梯度计算(gradient) [Fx,Fy,Fz,...] = gradient(F,h1,h2,...) [Fx,Fy,Fz,...] = gradient(F)
F是一个矩阵，F的维数表示自变量的个数n，返回n个大小相同的矩阵
v = -2:0.2:2;
[x,y] = meshgrid(v);
z = x .* exp(-x.^2 -y.^2);
[px,py] = gradient(z,.2,.2);
contour(v,v,z), hold on, quiver(v,v,px,py), hold off
第四章 数值积分与数值微分 1.确定下列求积公式中的特定参数，使其代数精度尽量高，并指明所构造出的求积公式所具有的代数精度： 101210121 12120 (1)()()(0)(); (2)()()(0)(); (3)()[(1)2()3()]/3; (4)()[(0)()]/2[(0)()]; h h h h h f x dx A f h A f A f h f x dx A f h A f A f h f x dx f f x f x f x dx h f f h ah f f h -----≈-++≈-++≈-++''≈++-?? ?? 解： 求解求积公式的代数精度时，应根据代数精度的定义，即求积公式对于次数不超过m 的多项式均能准确地成立，但对于m+1次多项式就不准确成立，进行验证性求解。 （1）若101(1) ()()(0)()h h f x dx A f h A f A f h --≈-++? 令()1f x =，则 1012h A A A -=++ 令()f x x =，则 110A h A h -=-+ 令2 ()f x x =，则 3 221123 h h A h A -=+ 从而解得 011431313A h A h A h -?=?? ? =?? ?=?? 令3 ()f x x =，则 3()0h h h h f x dx x dx --==? ? 101()(0)()0A f h A f A f h --++=
令4 ()f x x =，则 455 1012()5 2 ()(0)()3 h h h h f x dx x dx h A f h A f A f h h ---== -++=? ? 故此时， 101()()(0)()h h f x dx A f h A f A f h --≠-++? 故 101()()(0)()h h f x dx A f h A f A f h --≈-++? 具有3次代数精度。 （2）若 21012()()(0)()h h f x dx A f h A f A f h --≈-++? 令()1f x =，则 1014h A A A -=++ 令()f x x =，则 110A h A h -=-+ 令2 ()f x x =，则 3 2211163 h h A h A -=+ 从而解得 1143 8383A h A h A h -?=-?? ? =?? ?=?? 令3 ()f x x =，则 22322()0h h h h f x dx x dx --==? ? 101()(0)()0A f h A f A f h --++=
实验三 数值积分程序设计算法 1）实验目的 通过本次实验熟悉并掌握各种数值积分算法及如何在matlab 中通过设计程序实现这些算法，从而更好地解决实际中的问题。 2）实验题目 给出积分 dx x I ? -= 3 2 2 1 1 1.用Simpson 公式和N=8的复合Simpson 公式求积分的近似值. 2.用复合梯形公式、复合抛物线公式、龙贝格公式求定积分，要求绝对误差为 7 10*2 1-= ε，将计算结果与精确解做比较，并对计算结果进行分析。 3）实验原理与理论基础 Simpson 公式 )]()2 ( 4)([6 b f b a f a f a b S +++-= 复化梯形公式 将定积分? = b a dx x f I )(的积分区间],[b a 分隔为n 等分，各节点为 n j jh a x j ,,1,0, =+= n a b h -= 复合梯形(Trapz)公式为 ])()(2)([21 1 ∑-=++-= n j j n b f x f a f n a b T 如果将],[b a 分隔为2n 等分，而n a b h /)(-=不变， 则 )]()(2)(2)([41 2 111 2b f x f x f a f n a b T n j j n j j n +++-= ∑∑-=+-= 其中 h j a h x x j j )2 1(2 12 1+ +=+ =+ ，)]()(2)(2)([41 2 11 1 2b f x f x f a f n a b T n j j n j j n +++-= ∑∑-=+ -= ∑ -=-++-+ =1 )2) 12((22 1n j n n a b j a f n a b T n=1时，a b h -=，则)]()([2 1b f a f a b T +-= )0(0T = )2 1(2 2 112h a f a b T T + -+ =)1(0T = 若12-=k n ，记)1(0-=k T T n ， ,2,1=k 1 2 --= k a b h jh a x j +=1 2 --+=k a b j a h x x j j 2 12 1+ =+ k a b j a 2 ) 12(-++=，则可得如下递推公式
M a t l a b数值积分与数值微分 Matlab数值积分 1.一重数值积分的实现方法 变步长辛普森法、高斯-克朗罗德法、梯形积分法 1.1变步长辛普森法 Matlab提供了quad函数和quadl函数用于实现变步长 辛普森法求数值积分.调用格式为: [I,n]=Quad(@fname,a,b,tol,trace) [I,n]=Quadl(@fname,a,b,tol,trace) Fname是函数文件名，a,b分别为积分下限、积分上限； tol为精度控制，默认为1.0×10-6，trace控制是否展 开积分过程，若为0则不展开，非0则展开，默认不展开. 返回值I为积分数值；n为调用函数的次数. --------------------------------------------------------------------- 例如：求 ∫e0.5x sin(x+π )dx 3π 的值. 先建立函数文件 fesin.m function f=fesin(x) f=exp(-0.5*x).*sin(x+(pi/6));再调用quad函数
[I,n]=quad(@fesin,0,3*pi,1e-10) I= 0.9008 n= 365 --------------------------------------------------------------------- 例如：分别用quad函数和quadl函数求积分 ∫e0.5x sin(x+π 6 )dx 3π 的近似值，比较函数调用的次数. 先建立函数文件 fesin.m function f=fesin(x) f=exp(-0.5*x).*sin(x+(pi/6)); formatlong [I,n]=quadl(@fesin,0,3*pi,1e-10) I= n= 198 [I,n]=quad(@fesin,0,3*pi,1e-10) I= n= 365 --------------------------------------------------------------------- 可以发现quadl函数调用原函数的次数比quad少，并 且比quad函数求得的数值解更精确. 1.2高斯-克朗罗德法
专题六数值微积分与方程求解6.1 数值微分与数值积分 ?数值微分 ?数值积分
1．数值微分 （1）数值差分与差商 微积分中，任意函数f(x)在x 0点的导数是通过极限定义的： h x f h x f x f h )()(lim )('0 0-+=→h h x f x f x f h ) ()(lim )('0 00 0--=→h h x f h x f x f h ) 2/()2/(lim )('0 --+=→
) ()()(000 x f h x f x f -+=?) ()()(0 h x f x f x f --=?) 2/()2/()(0 h x f h x f x f --+=δ如果去掉极限定义中h 趋向于0的极限过程，得到函数在x 0点处以h （h>0）为步长的向前差分、向后差分和中心差分公式： 向前差分： 向后差分： 中心差分：
函数f(x)在点x 0的微分接近于函数在该点的差分，而f 在点x 的导数接近于函数在该点的差商。 h x f h x f x f ) ()(≈ )('0 00 -+h h x f x f x f ) ()(≈ )('0 00 --h h x f h x f x f ) 2/()2/(≈ )('0 --+向前差商： 向后差商： 中心差商： 当步长h 充分小时，得到函数在x 0点处以h （h>0）为步长的向前差商、 向后差商和中心差商公式：
（2）数值微分的实现 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，按行计算差分。 注意：diff函数计算的是向量元素间的差分，故差分向量元素的个数比原向量少了一个。同样，对于矩阵来说，差分后的矩阵比原矩阵少了一行或一列。 另外，计算差分之后，可以用f(x)在某点处的差商作为其导数的近似值。
数值积分与 数值微分 习题课
一、已知012113,,424x x x ===，给出以这 3个点为求积节 点在[]0.1上的插值型求积公式 解：过这3个点的插值多项式基函数为 ()()()()()()()()()()()()()()()()1202010202121012012220211 20,0,1,2 k k x x x x l x x x x x x x x x l x x x x x x x x x l x x x x x A l x dx k --= ----= ----= --==?
()()()()()()()()()()()()111200001021102100101210120202113224111334244131441113324241142x x x x x x A dx dx x x x x x x x x x x A dx dx x x x x x x x x x x A dx x x x x ????-- ???--????=== --????-- ??? ???? ????-- ???--????===- --????-- ??? ???? ????-- ??--???==--?????102313134442dx ??= ????-- ??? ???? ? 故所求的插值型求积公式为 ()1 211 123343234f x dx f f f ??????≈- + ? ? ??????? ?
二、确定求积公式 ( )( )(1 1158059f x dx f f f -? ?≈++?? ? 的代数精度，它是Gauss 公式吗？ 证明：求积公式中系数与节点全部给定，直接检验 依次取()23451,,,,,f x x x x x x =，有 [ ](1 1111 215181519 1058059dx xdx --==?+?+???==?+?+?? ???
第4章附录 4.2.2 复化求积分 例题4.2.5计算程序 //simp.c// # include
subroutine trap(a,b,f,eps,t,n) real(8) a,b,f,t,fa,fb,h,t0,s,x fa=f(a); fb=f(b) n=1; h=b-a t0=h*(fa+fb)/2.0 5 s=0.0 do 10 k=0,n-1 x=a+(k+0.5)*h s=s+f(x) 10 continue t=(t0+h*s)/2.0 if (abs(t-t0).ge.eps) then t0=t n=n+n h=h/2.0 goto 5 end if return end %%% demo_aTrapInt.m %%% function demo_aTrapInt clc;clear; format long; [T nsub] = aTrapInt(@f01,0,1,0.000001) end function [T nsub]= aTrapInt(f,a,b,eps) tol = 1; nsub = 1; inall = 0; T = 0.5*(b-a)*(f(a)+f(b)); while tol > eps T0 = T; nsub = 2*nsub; n = nsub + 1; % total number of nodes h = (b-a)/nsub; % stepsize x = a:h:b; % divide the interval inall = inall + sum(f(x(2:2:n-1))); T = 0.5*h * (f(a)+2*inall+f(b)); tol = abs(T-T0); end end
第五章 数值微分与数值积分 一．分别用向前差商，向后差商和中心差商公式计算()f x =2x =的导数的近似值。其中，步长0.1h =。 【详解】 00()()(20.1)(2)=0.349 2410.10.1 f x h f x f f h +?+?===向前差商 00()()(2)(20.1)=0.358 0870.10.1 f x f x h f f h ????===向后差商 00()()(20.1)(20.1)= 0.353 664220.10.2f x h f x h f f h +??+??===×中心差商 二．已知数据 x 2.5 2.55 2.60 2.65 2.70 ()f x 1.58114 1.59687 2 1.62788 1.64317 求( 2.50),(2.60),(2.70)f f f ′′′的近似值。 【详解】 0.05h =，按照三点公式 3(2.50)4(2.55)(2.60)3 1.581144 1.59687 1.61245(2.50)0.316 10020.050.1 f f f f ?+??×+×?′≈==×(2.65)(2.55)1.627881.59687(2.60)0.310 10020.050.1 f f f ??′≈==× (2.60)4(2.65)3(2.70)241.6278831.64317(2.70) 4.179 90020.050.1 f f f f ?+?×+×′≈==× 三．已知如下数据 x 3 4 5 6 7 8 ()f x 2.937 6 6.963 213.600 0 23.500 8 37.318 4 55.705 6
数值分析实验报告四 数值积分与数值微分实验（2学时） 一 实验目的 1．掌握复化的梯形公式、Simpson 公式等牛顿-柯特斯公式计算积分。 2. 掌握数值微分的计算方法。 二 实验内容 1. 用复化梯形公式计算积分。 ?9 0dx x M=8 2. 用复化Simpson 公式计算积分。 ? 90dx x M=8 3. 给定下列表格值 利用四点式（n=3）求)50()50('''f f 和的值。 三 实验步骤（算法）与结果 1复化梯形公式 用C 语言编程如下： #include
float y; y=sqrt(x); return y; } void main() { int i,m; float a,b,h,r; printf("输入等分数m:" ); scanf("%d",&m); printf("输入区间左端点a的值:"); scanf("%f",&a); printf("输入区间右端点b的值:"); scanf("%f",&b); float x[m+1]; h=(b-a)/m; for(i=0;i<=m;i++) x[i]=a+i*h; r=0; for(i=0;i<=m;i++) {if(i==0) r=r+h*0.5*f(x[i]); if(i>0&&i if(i==m) r=r+0.5*h*f(x[i]); } printf("输出区间[%3.1f %3.1f]的积分值:%f\n",a,b,r); } 求解结果如下： 输入等分数m:8 输入区间左端点a的值:0 输入区间右端点b的值:9 输出区间[0.0 9.0]的积分值:17.769514 2复化Simpson公式 用C语言编程如下： #include 数值积分与微分 摘要 本文首先列举了一些常用的数值求积方法，一是插值型求积公式，以N e w t o n C o t e s -公式为代表，并分析了复合型的Newton Cotes -公式；另一个是Gauss Ledendre -求积公式，并给出几个常用的Gauss Ledendre -求积公式。其次，本文对数值微分方法进行分析，主要是差分型数值微分和插值型数值微分，都给出了几种常用的微分方法。然后，本文比较了数值积分与微分的关系，发现数值积分与微分都与插值或拟合密不可分。 本文在每个方法时都分析了误差余项，并且在最后都给出了MATLAB 的调用程序。 关键词：插值型积分Gauss Ledendre -差分数值微分插值型数值微分 MATLAB 一、常用的积分方法 计算积分时，根据Newton Leibniz -公式， ()()()b a f x dx F b F a =-? 但如果碰到以下几种情况： 1）被积函数以一组数据形式表示； 2）被积函数过于特殊或者原函数无法用初等函数表示 3）原函数十分复杂难以计算 这些现象中，Newton Leibniz -公式很难发挥作用，只能建立积分的近似计算方法，数值积分是常用的近似计算的方法。 1.1 插值型积分公式 积分中的一个常用方法是利用插值多项式来构造数值求积公式，具体的步骤如下： 在积分区间上[,]a b 上取一组节点：01201,,,,()n n x x x x a x x x b ≤<<≤ 。已知()k x f 的函数值，作()x f 的n 次插值多项式，则 (1) ()10()()()()() (1)!n n x n n k k n k f f L x R x f x l x w x n ++==+=++∑ 其中，()k l x 为n 次插值基函数，则得 (1)+10()(()())1 =[()]()[()](1)!b b n n a a n b b n k k n a a k f x dx L x R x dx l x dx f x f x w x dx n ξ+==+++? ?∑??（） 公式写成一般形式： ()()[]n b k k n a k f x dx A f x R f ==+∑? 其中， 01100110 ()()()() ()()()()()b b k k k k a a k k k k k k x x x x x x x x A l x dx dx x x x x x x x x -+-+----==----?? (1)+11 [][()](1)!b n n n a R f f x w x dx n ξ+=+?（） 显然，当被积函数f 为次数小于等于n 的多项式时，其相应的插值型求积公式为准确公式，即： ()() n b k k a k f x dx A f x ==∑? 1.1.1 求积公式的代数精度 定义：求积公式对于任何次数不大于m 的代数多项式()f x 均精确成立，而对于 1()m f x x +=不精确成立，则称求积公式具有m 次代数精度。 定理：含有1n +个节点(0,1,,)k x k n = 的插值型求积公式的代数精度至少为n 。 第七章数值积分与数值微分 积分问题最早来自于几何形体的面积、体积计算，也是经典力学中的重要问题（例如计算物体的重心位置）. 在现实应用中，很多积分的结果并不能写成解析表达式，因此需要通过数值方法来计算. 数值微分是利用一些离散点上的函数值近似计算某一点处的函数导数，它针对表达式未知的函数. 本章介绍一元函数积分（一重积分）和微分的各种数值算法，它们也是数值求解积分方程、微分方程的基础. 7.1数值积分概论 7.1.1基本思想 考虑如下定积分的计算： I(f)≡∫f(x)dx b a ,(7.1) 其中函数f: ?→?，首先应想到的是微积分中学习过的牛顿-莱布尼兹（Newton-Leibniz）公式： ∫f(x)dx b a =F(b)?F(a) , 其中F′(x)=f(x),即F(x)为f(x)的原函数. 但是，诸如e x2,sinx x ,sinx2等表达式很简单的函数却找不到用初等函数表示的原函数，因此必须研究数值方法来近似计算积分. 另一方面，某些函数的原函数虽然可以解析表示，但其推导、计算非常复杂，此时也需要使用数值积分方法. 一般考虑连续的、或在区间[a,b]上可积①的函数f(x)，则根据积分的定义有： lim n→∞,?→0∑(x i+1?x i)f(ξi) n i=0 =I(f) , (7.2) 其中a=x0 第5章数值积分与数值微分方法关于定积分计算，已经有较多方法，如公式法、分步积分法等，但实际问题中，经常出现不能用通常这些积分方法计算的定积分问题。怎样把这些通常方法失效的定积分在一定精度下快速计算出来，特别是通过计算机编程计算出来就是本章研究的内容。 此外，怎样根据函数在若干个点处的函数值去求该函数的导数近似值也是本章介绍的内容。 本章涉及的方法有Newton-Cotes求积公式、Gauss求积公式、复化求积公式、Romberg求积公式和数值微分。 5.1 引例 人造地球卫星轨道可视为平面上的椭圆。我国的第一颗人造地球卫星近地点距离地球表面439km ，远地点距地球表面2384km ，地球半径为6371km ，求该卫星的轨道长度。 本问题可用椭圆参数方程 cos ,,0sin x a t a b y b t π=?≤≤>?=? (0t 2) 来描述人造地球卫星的轨道，式中a, b 分别为椭圆的长短轴，该轨道的长度L 就是如下参数方程弧长积分 但这个积分是椭圆积分，不能用解析方法计算。 5.2问题的描述与基本概念 要想用计算机来计算()b a f x dx ?，应对其做离散化处 理。注意到定积分是如下和式的极限 1 ()lim ()n b i i a i f x dx f x λξ→==?∑? 要离散化，做 1） 去掉极限号lim 2） 将i ξ取为具体的i x 值 3） 为减少离散化带来的误差，将i x ?用待定系数i A 代替 于是就得到 定义 5.1 若存在实数1212,,,;,, ,,n n x x x A A A 且任 取()[,],f x C a b ∈都有 1 ()()n b i i a i f x dx A f x =≈∑? （5.1） 2.3 数值积分 2.3.1 一元函数的数值积分 函数1 quad 、quadl 、quad8 功能 数值定积分，自适应Simpleson 积分法。 格式 q = quad(fun,a,b) %近似地从a 到b 计算函数fun 的数值积分，误差为10-6。 若给fun 输入向量x ，应返回向量y ，即fun 是一单值函数。 q = quad(fun,a,b,tol) %用指定的绝对误差tol 代替缺省误差。tol 越大，函数计 算的次数越少，速度越快，但结果精度变小。 q = quad(fun,a,b,tol,trace,p1,p2,…) %将可选参数p1,p2,…等传递给函数 fun(x,p1,p2,…)，再作数值积分。若tol=[]或trace=[]，则用缺省值进行计算。 [q,n] = quad(fun,a,b,…) %同时返回函数计算的次数n … = quadl(fun,a,b,…) %用高精度进行计算，效率可能比quad 更好。 … = quad8(fun,a,b,…) %该命令是将废弃的命令，用quadl 代替。 例2-40 >>fun = inline(‘3*x.^2./(x.^3-2*x.^2+3)’); equivalent to: function y=funn(x) y=3*x.^2./(x.^3-2*x.^2+3); >>Q1 = quad(fun,0,2) >>Q2 = quadl(fun,0,2) 计算结果为： Q1 = 3.7224 Q2 = 3.7224 补充：复化simpson 积分法程序 程序名称 Simpson.m 调用格式 I=Simpson('f_name',a,b,n) 程序功能 用复化Simpson 公式求定积分值 输入变量 f_name 为用户自己编写给定函数()y f x 的M 函数而命名的程序文件名 a 为积分下限 b 为积分上限 n 为积分区间[,]a b 划分成小区间的等份数 输出变量 I 为定积分值 程序 function I=simpson(f_name,a,b,n) h=(b-a)/n; x=a+(0:n)*h; f=feval(f_name,x); 解：卫星轨道的示意图如右上图所示，,a b 分别是椭圆轨道的长半轴和短半轴，地球位于椭圆的一个焦点处，焦距为c ，地球半径为r ，近地点和远地点与地球表面的距离分别是1s 和2s . 由图中可知，上述数据存在如下关系： 12122,a s s r c a s r =++=-- 由椭圆性质 b =，将12,,s s r 的数据代入以上各式可得7782.5a km =，7721.5b km =. 椭圆的参数方程为： c o s ,s i n x a t y b t == ， (02)t π≤< 根据计算参数方程弧长的公式，椭圆长度可表为如下积分： /2 22221/20 4(sin cos )L a t b t dt π=+? 由于该积分无法求得解析解，下面我们编写MATLAB 程序对其进行数值求解。 s1=439;s2=2384;r=6371; a=(s1+s2)/2+r a = 7.7825e+003 >> c=a-s1-r; >> b=sqrt(a^2-c^2) b = 7.7215e+003 y=inline('sqrt(7782.5^2*sin(t).^2+7721.5^2*cos(t).^2)'); %建立积分内联函数 >> t=0:pi/10:pi/2; y1=y(t); format long >> L1=4*trapz(t,y1) %梯形积分 L1 = 4.870744099902405e+004 >> L2=4*quad(y,0,pi/2,1e-6) %辛普森积分 L2 = 4.870744099903280e+004 求解结果显示：两种方法计算求得的积分结果相当接近，轨道长度约为：4 4.8710km ?. 第4章数值积分与数值微分 计算定积分是科学技术和工程计算中经常遇到的问题，然而绝大多数求定积分的问题都不能通过被积函数的原函数来求解，而要用数值的方法来解决. 4．1问题的提出神舟十号载人飞船于2013年6月11日17时38分在酒泉卫星发射中心成功发射，6月13日13时18分与天宫一号成功对接，在轨飞行15天，其中12天与天宫一号组成组合体在太空中飞行，6月26日8时7分，神舟十号返回舱成功返回地面. 神舟十号载人飞船发射的初始轨道为近地点约200公里、远地点约330公里的椭圆轨道，对接轨道为距地约343公里的近圆轨道，飞行速度约为每秒7.9公里.试计算神舟十号载人飞船在椭圆轨道飞行一圈的公里数，进而可以计算在轨飞行的公里数. 这里主要是计算神舟十号的椭圆轨道的周长.如图4.1所示， 图4.1神舟十号的椭圆轨道 已知地球半径=r 6371km，近地点高度为=1s 200km、远地点高度为=2s 330km，则椭圆长半轴=++=2/)2(21s s r a 6636km ，半焦距=-=2/)(12s s c 65km ，由椭圆参数方程t b y t a x sin ,cos ==，其中22c a b -=，知椭圆周长为 dt t b t a L ?+=20 2222cos sin 4π dt t c a ?-=20 222cos 4π dt t a c a ?-=2 0222cos 14π 这是一个定积分，只要求出它的值就行了. 一般的，对于定积分 ()b a I f x dx =?如果被积函数()f x 在区间[,]a b 上连续，且()f x 的原函数为()F x ，则由牛顿—莱布尼兹（Newton-Leibnitz）公式，有 ()()() b a f x dx F b F a =-?似乎问题已经解决，其实不然. （1）有很多被积函数找不到用初等函数表示的原函数F(x)，例如 等等，表面看它们并不复杂，但却无法求得用初等函 数表示的原函数F(x).这里的椭圆积分dt t a c ?-2 0222cos 1π 也是这样，（2）有的积分即使能找到用初等函数表示的原函数F(x)，但原函数非常复杂，用牛顿—莱布尼兹公式，计算也很困难，例如 24211211ln arc (21)arc (21)1422122 x x dx tg x tg x C x x x ++??=+++-+??+-+?（3）有的被积函数()f x 是由测量或数值计算给出的数据表，是列表函数，也无法用牛顿—莱布尼兹公式计算. 对于上述这些情况，都要求建立定积分的近似计算方法——数值积分法.数值积分与微分方法
第7章 数值积分与数值微分
数值积分与数值微分
数值积分与微分方程
数值积分与数值微分 编程计算
数值积分与数值微分
相关主题
文本预览