资源描述:
《刘卫国全套配套课件MATLAB程序设计与应用第3版第7章 MATLAB数值微分与积分__源程序.doc》由会员上传分享,免费在线阅读,更多相关内容在教育资源-天天文库。
1、第7章MATLAB数值微分与积分例7-1设x由[0,2π]间均匀分布的6个点组成,求sinx的1~3阶差分。>>X=linspace(0,2*pi,6);>>Y=sin(X)>>DY=diff(Y)%计算Y的一阶差分>>D2Y=diff(Y,2)%计算Y的二阶差分,也可用命令diff(DY)计算>>D3Y=diff(Y,3)%计算Y的三阶差分,也可用diff(D2Y)或diff(DY,2)例7-2用不同的方法求函数f(x)的数值导数,并在同一个坐标系中做出f'(x)的图像。f=@(x)sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2;g=@(x)(
2、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求导数dpdpx=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,'-')%作图例7-3求。先建立一个函数文件fex.m。functionf=fex(
3、x)f=exp(-x.^2);接下来调用数值积分函数quad求定积分,命令如下:>>[I,n]=quad(@fex,0,1)也可不建立关于被积函数的函数文件,而使用匿名函数(或内联函数)求解,命令如下:>>f=@(x)exp(-x.^2);%用匿名函数f(x)定义被积函数>>[I,n]=quad(f,0,1)%注意函数句柄f前面不加@号例7-4分别用quad函数和quadl函数求积分的近似值,并在相同的积分精度下,比较函数的调用次数。>>formatlong4>>f=@(x)4./(1+x.^2);>>[I,n]=quad(f,0,1,1e-8)%调用函数quad求定积分>>[I,
4、n]=quadl(f,0,1,1e-8)%调用函数quadl求定积分>>(atan(1)-atan(0))*4%理论值>>format例7-5求。先建立被积函数文件fe.m。functionf=fe(x)f=1./(x.*sqrt(1-log(x).^2));再调用数值积分函数integral求定积分,命令如下:>>I=integral(@fe,1,exp(1))例7-6求。建立被积函数文件fsx.m。functionf=fsx(x)f=sin(1./x)./x.^2;调用函数quadgk求定积分,命令如下:>>I=quadgk(@fsx,2/pi,+Inf)例7-7用trapz函
5、数计算定积分。>>formatlong>>x=0:0.001:1;>>y=4./(1+x.^2);%生成函数向量>>trapz(x,y)>>format5.累计梯形积分在MATLAB中,提供了对数据积分逐步累计的函数cumtrapz。该函数调用格式如下。Z=cumtrapz(Y)Z=cumtrapz(X,Y)函数其他参数的含义和用法与trapz函数的相同。例如:>>S=cumtrapz([1:5;2:6]')例7-8计算二重定积分。4建立一个函数文件fxy.m。functionf=fxy(x,y)globalki;ki=ki+1;%ki用于统计被积函数的调用次数f=exp(-x.^
6、2/2).*sin(x.^2+y);调用函数求解,命令如下:>>globalki;>>ki=0;>>I=integral2(@fxy,-2,2,-1,1)%调用integral2函数求解>>ki=0;>>I=quad2d(@fxy,-2,2,-1,1)%调用quad2d函数求解>>ki>>ki=0;>>I=dblquad(@fxy,-2,2,-1,1)%调用dblquad函数求解例7-9计算三重定积分>>fxyz=@(x,y,z)4*x.*z.*exp(-z.*z.*y-x.*x);%定义被积函数>>integral3(fxyz,0,pi,0,pi,0,1)>>triplequad
7、(fxyz,0,pi,0,pi,0,1,1e-7)例7-10给定数学函数x(t)=12sin(2π×10t+π/4)+5cos(2π×40t)取N=128,试对t从0~1s采样,用fft函数作快速傅里叶变换,绘制相应的振幅—频率图。程序如下:N=128;%采样点数T=1;%采样时间终点t=linspace(0,T,N);%给出N个采样时间x=12*sin(2*pi*10*t+pi/4)+5*cos(2*pi*40*t);%求各采样点样本值xdt=t(2)-t(1);%