Matlab数值积分与数值微分

Matlab数值积分与数值微分
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π
的近似值，⽐较函数调⽤的次数.
先建⽴函数⽂件
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⾼斯-克朗罗德法
Matlab提供了⾃适应⾼斯-克朗罗德法的quadgk函数来求震荡函数的定积分，函数的调⽤格式为：
[I,err]=quadgk(@fname,a,b)
Err返回近似误差范围，其他参数的意义与quad函数相同，积分上下限可以是-Inf或Inf，也可以是复数，若为复数则在复平⾯上求积分.
--------------------------------------------------------------------- 例如：求积分
∫
xsinx
1+cos2x
dx π
的数值.
先编写被积函数的m⽂件
fsx.m
function f=fsx(x)
f=x.*sin(x)./(1+cos(x).^2);
再调⽤quadgk函数
I=quadgk(@fsx,0,pi)
I=
2.4674
--------------------------------------------------------------------- 例如：求积分
∫
xsinx
dx +∞
∞
的值.
先编写被积函数的m⽂件
fsx.m
function f=fsx(x)
f=x.*sin(x)./(1+cos(x).^2); 再调⽤quadgk函数
I=quadgk(@fsx,-Inf,Inf)
I=
-9.0671e+017
---------------------------------------------------------------------
1.3梯形积分法
对于⼀些不知道函数关系的函数问题，只有实验测得的⼀
组组样本点和样本值，由表格定义的函数关系求定积分问
题⽤梯形积分法，其函数是trapz函数，调⽤格式为：
I=Traps(X,Y)
X,Y为等长的两组向量，对应着函数关系Y=f(X) X=(x1,x2,…,x n)(x1
分区间是[x1,x n]
--------------------------------------------------------------------- 例如：已知某次物理实验测得如下表所⽰的两组样本点.
现已知变量x和变量y满⾜⼀定的函数关系，但此关系
未知，设y=f(x)，求积分
13.39
∫f(x)dx
1.38
的数值.
X=[1.38,1.56,2.21,3.97,5.51,7.79,9.19,11
.12,13.39];
Y=[3.35,3.96,5.12,8.98,11.46,17.63,24.41
,29.83,32.21]; I=trapz(X,Y) I=
217.1033
---------------------------------------------------------------------
例如：⽤梯形积分法求积分：
∫e ?x dx 2.5
1
的数值.
x=1:0.01:2.5; y=exp(-x); I=trapz(x,y) I= 0.2858
---------------------------------------------------------------------
2. 多重数值积分的实现
重积分的积分函数⼀般是⼆元函数f(x,y)或三元函数f(x,y,z)；形如：
∫∫f (x,y )dxdy b
a d
c
∫∫∫f(x,y,z)dxdydz b a d c
f e
Matlab 中有dblquad 函数和triplequad 函数来对上述两个积分实现.调⽤格式为： I=dblquad(@fun,a,b,c,d,tol)
I=triplequad(@fun,a,b,c,d,e,f,tol)
Fun 为被积函数，[a,b]为x 的积分区间；[c,d]为y 的积分区间；[e,f]为z 的积分区间.
Dblquad 函数和triplequad 函数不允许返回调⽤的次数，如果需要知道函数调⽤的次数，则在定义被积函数的m ⽂件中增加⼀个计数变量，统计出被积函数被调⽤的次数.
---------------------------------------------------------------------
例如:计算⼆重积分
I =∫∫√dxdy π2
π2π2
π2
的值.
先编写函数⽂件fxy.m
function f=fxy(x,y) global k; k=k+1;
f=sqrt(x.^2+y.^2);
再调⽤函数dblquad
globalk; k=0;
I=dblquad(@fxy,-pi/2,pi/2,-pi/2,pi/2,1.0e-10) I= 11.8629 k k= 37656
---------------------------------------------------------------------
例如：求三重积分
∫∫∫4xze ?z
2y?x 2
dxdydz π
π
1
的值.
编写函数⽂件fxyz1.m
function f=fxyz1(x,y,z)
global j;
j=j+1;
f=4*x.*z.*exp(-z.*z.*y-x.*x);
调⽤triplequad函数
edit
globalj;
j=0;
I=triplequad(@fxyz1,0,pi,0,pi,0,1,1.0e-10)
I=
1.7328
j
j=
1340978
---------------------------------------------------------------------Matlab数值微分
1.数值微分与差商
导数的三种极限定义
f′(x)=lim
n→0f(x+h)?f(x)
h
f′(x)=lim
n→0
f(x)?f(x?h)
f′(x)=lim
n→0f(x+h2)?f(x?h2)
h
上述公式中假设h>0，引进记号：
f(x)=f(x+h)f(x)
f(x)= f(x)f(xh)
δf(x)= f(x+h
)?f(x?
h
)
称上述?f(x)、?f(x)、δf(x)为函数在x点处以h(h>0)为步长的向前差分、向后差分、中⼼差分，当步长h⾜够⼩时，有：
f′(x)≈?f(x) h
f′(x)≈
f(x) f′(x)≈
δf(x)
f(x) h 、?f(x)
h
、δf(x)
h
也分别被称为函数在x点处以h(h>0)为
步长的向前差商、向后差商、中⼼差商.当h⾜够⼩时，函数f(x)在x点处的导数接近于在该点的任意⼀种差商，微分接近于在该点的任意⼀种差分.
2.函数导数的求法
2.1⽤多项式或样条函数g(x)对函数f(x)进⾏逼近(插
值或拟合)，然后⽤逼近函数g(x)在点x处的导数作
为f(x)在该点处的导数.
2.2⽤f(x)在点x处的差商作为其导数.
3.数值微分的实现⽅法
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按⾏计算差分.
--------------------------------------------------------------------- 例如：⽣成6阶范德蒙德矩阵，然后分别按⾏、按列计算
⼆阶向前差分
A=vander(1:6)
A=
111111
32168421
2438127931
1024256641641
31256251252551
777612962163661
D2A1=diff(A,2,1)
D2A1=
1805012200
57011018200
132019424200
255030230200
D2A2=diff(A,2,2)
D2A2=
0000
8421
10836124
576144369
20004008016
540090015025
--------------------------------------------------------------------- 例如：设f(x)=√x3+2x2?x+12+√(x+5)
6+5x+2求函数f(x)的数值导数，并在同⼀坐标系中作出f’(x)
的图像.
已知函数f(x)的导函数如下：
f′(x)=
3x2+4x?1
2√x3+2x2?x+12
+
1
6√()5
6
+5
编辑函数⽂件fun7.m和fun8.m functionf=fun7(x)
f=sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2;
functionf=fun8(x)
f=(3*x.^2+4*x-1)/2./sqrt(x.^3+2*x.^2-x+12)+1/6./(x+5).^(5/6)+5 ;
x=-3:0.01:3;
p=polyfit(x,fun7(x),5);⽤5次多项式拟合曲线
dp=polyder(p);对拟合多项式进⾏求导
dpx=polyval(dp,x);对dp在假设点的求函数值
dx=diff(fun7([x,3.01]))/0.01;直接对dx求数值导数
gx=fun8(x);求函数f的函数在假设点的导数
plot(x,dpx,x,dx,'.',x,gx,'-')
可以发现，最后得到的三条曲线基本重合.
--------------------------------------------------------------------- 练习：
A.⽤⾼斯-克朗罗德法求积分
∫
dx
1+x2 +∞
∞
的值并讨论计算⽅法的精确度.（该积分值为π）function f=fun9(x)
f=1./(1+x.^2);
formatlong
[I,err]=quadgk(@fun9,-Inf,Inf)
I=
err=
B.设函数
f(x)=
sin x
⽤不同的办法求该函数的数值导数，并在同⼀坐标系中作出f′(x)的图像.已知
f′(x)=x cos x+cos x cos2x?sin x+2sin x sin2x
()2
function f=fun10(x)
f=sin(x)./(x+cos(2*x));
function f=fun11(x)
f=(x.*cos(x)+cos(x).*cos(2*x)-sin(x)-2*sin(x).*sin(2*x))/(x+cos(2 *x)).^2; x=-3:0.01:3;
p=polyfit(x,fun10(x),5);
dp=polyder(p);
dpx=polyval(dp,x);
dx=diff(fun10([x,3.01]))/0.01;
gx=fun11(x);
plot(x,dpx,'r:',x,dx,'.g',x,gx,'-k')。