Matlab数值积分与数值微分.docx
- 文档编号:26816914
- 上传时间:2023-06-23
- 格式:DOCX
- 页数:12
- 大小:67.69KB
Matlab数值积分与数值微分.docx
《Matlab数值积分与数值微分.docx》由会员分享,可在线阅读,更多相关《Matlab数值积分与数值微分.docx(12页珍藏版)》请在冰豆网上搜索。
Matlab数值积分与数值微分
Matlab数值积分与数值微分
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为调用函数的次数.
---------------------------------------------------------------------
例如:
求
的值.
先建立函数文件
fesin.m
functionf=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函数求积分
的近似值,比较函数调用的次数.
先建立函数文件
fesin.m
functionf=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,也可以是复数,若为复数则在复平面上求积分.
---------------------------------------------------------------------
例如:
求积分
的数值.
先编写被积函数的m文件
fsx.m
functionf=fsx(x)
f=x.*sin(x)./(1+cos(x).^2);
再调用quadgk函数
I=quadgk(@fsx,0,pi)
I=
2.4674
---------------------------------------------------------------------
例如:
求积分
的值.
先编写被积函数的m文件
fsx.m
functionf=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,…,xn)(x1 --------------------------------------------------------------------- 例如: 已知某次物理实验测得如下表所示的两组样本点. 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 现已知变量x和变量y满足一定的函数关系,但此关系未知,设y=f(x),求积分 的数值. 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 --------------------------------------------------------------------- 例如: 用梯形积分法求积分: 的数值. x=1: 0.01: 2.5; y=exp(-x); I=trapz(x,y) I= 0.2858 --------------------------------------------------------------------- 2.多重数值积分的实现 重积分的积分函数一般是二元函数f(x,y)或三元函数f(x,y,z);形如: 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文件中增加一个计数变量,统计出被积函数被调用的次数. --------------------------------------------------------------------- 例如: 计算二重积分 的值. 先编写函数文件fxy.m functionf=fxy(x,y) globalk; 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 --------------------------------------------------------------------- 例如: 求三重积分 的值. 编写函数文件fxyz1.m functionf=fxyz1(x,y,z) globalj; 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.数值微分与差商 导数的三种极限定义 上述公式中假设h>0,引进记号: 称上述 为函数在x点处以h(h>0)为步长的向前差分、向后差分、中心差分,当步长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 132********0 255030230200 D2A2=diff(A,2,2) D2A2= 0000 8421 10836124 576144369 20004008016 540090015025 --------------------------------------------------------------------- 例如: 设 求函数f(x)的数值导数,并在同一坐标系中作出f’(x)的图像. 已知函数f(x)的导函数如下: 编辑函数文件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.用高斯-克朗罗德法求积分 的值并讨论计算方法的精确度.(该积分值为π) functionf=fun9(x) f=1./(1+x.^2); formatlong [I,err]=quadgk(@fun9,-Inf,Inf) I= err= B.设函数 用不同的办法求该函数的数值导数,并在同一坐标系中作出 的图像. 已知 functionf=fun10(x) f=sin(x)./(x+cos(2*x)); functionf=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')
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Matlab 数值 积分 微分
![提示](https://static.bdocx.com/images/bang_tan.gif)