第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)
(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,'-')
例7-3 求I?-xe? 0dx。 12%作图
先建立一个函数文件fex.m。 function f=fex(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函数求积分的近似值，并在相同的积分精度下，比较函数的调用次数。
>> format long
4
>> f=@(x) 4./(1+x.^2);
>> [I,n]=quad(f,0,1,1e-8) %调用函数quad求定积分 >> [I,n]=quadl(f,0,1,1e-8) %调用函数quadl求定积分 >> (atan(1)-atan(0))*4 %理论值 >> format
例7-5 求I?e?x111?lnx2dx。
先建立被积函数文件fe.m。 function f=fe(x)
f=1./(x.*sqrt(1-log(x).^2));
再调用数值积分函数integral求定积分，命令如下： >> I=integral(@fe,1,exp(1))
例7-6 求
?? ??1 2x2sin1dx。 x建立被积函数文件fsx.m。 function f=fsx(x) f=sin(1./x)./x.^2;
调用函数quadgk求定积分，命令如下： >> I=quadgk(@fsx,2/pi,+Inf)
例7-7 用trapz函数计算定积分。 >> format long >> x=0:0.001:1;
>> y=4./(1+x.^2); %生成函数向量 >> trapz(x,y) >> format 5．累计梯形积分
在MATLAB中，提供了对数据积分逐步累计的函数cumtrapz。该函数调用格式如下。 Z=cumtrapz(Y) Z=cumtrapz(X,Y)
函数其他参数的含义和用法与trapz函数的相同。例如： >> S=cumtrapz([1:5;2:6]')
例7-8 计算二重定积分。
4
建立一个函数文件fxy.m。 function f=fxy(x,y) global ki;
ki=ki+1; %ki用于统计被积函数的调用次数 f=exp(-x.^2/2).*sin(x.^2+y); 调用函数求解，命令如下： >> global ki; >> 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 计算三重定积分
???001??04xze?z2y?x2dxdydz
>> fxyz=@(x,y,z) 4*x.*z.*exp(-z.*z.*y-x.*x); %定义被积函数 >> integral3(fxyz,0,pi,0,pi,0,1) >> triplequad(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); % 求各采样点样本值x dt=t(2)-t(1); % 采样周期 f=1/dt; % 采样频率(Hz)
X=fft(x); % 计算x的快速傅里叶变换X F=X(1:N/2+1); % F(k)=X(k)(k=1:N/2+1) f=f*(0:N/2)/N; % 使频率轴f从零开始 plot(f,abs(F),'-*') % 绘制振幅—频率图
4
xlabel('Frequency'); ylabel('|F(k)|')
求X的快速傅里叶逆变换，并与原函数进行比较，命令如下： >> ix=real(ifft(X)); %求逆变换，结果只取实部 >> plot(t,x,t,ix,':') %逆变换结果和原函数的曲线 >> norm(x-ix) %逆变换结果和原函数之间的距离
4
12454-915222-MATLAB程序设计与应用-第7章 MATLAB数值微分与积分__源程序
