1、第7章数值微积分与常微分方程求解【本章学习目标】掌握数值微分的实现方法。掌握数值积分基本原理和实现方法。了解离散傅里叶变换的原理及实现方法。掌握常微分方程的数值求解方法。7.1 数值微分称f(x)、f(x)及f(x)分别为函数在x点处以h（h0）为步长的向前差分、向后差分和中心差分。称f(x)/h、f(x)/h及f(x)/h分别为函数在x点处以h（h0）为步长的向前差商、向后差商和中心差商。7.1.1 数值差分与差商c1 1多项式求导法多项式求导法用多项式或样条函数g(x)对f(x)进行逼近（插值或拟合），然后用逼近函数g(x)在点x处的导数作为f(x)在点x处的导数。2用diff函数计算差分
2、用f(x)在点x处的某种差商作为其导数。在MATLAB中，计算向前差分的函数diff，调用格式如下。DX=diff(X)：计算向量X的向前差分。DX=diff(X,n)：计算向量X的n阶向前差分。DX=diff(X,n,dim)：计算矩阵X的n阶差分，dim=1时（默认状态），按列计算差分；dim=2，按行计算差分。7.1.2 数值微分的实现7.1 数值微分【例7.17.1】设f(x)=sinx，用不同的方法求函数f(x)的数值导数，并在同一个坐标系中绘制f(x)的三种方法所得导数曲线。x=0:pi/24:pi;%用5次多项式p拟合f(x),并对拟合多项式p求导数dp在假设点的函数值p=pol
3、yfit(x,sin(x),5);dp=polyder(p);dpx=polyval(dp,x);%直接对sin(x)求数值导数dx=diff(sin(x,pi+pi/24)/(pi/24);%求函数f的导函数g在假设点的导数gx=cos(x);plot(x,dpx,b-,x,dx,ko,x,gx,r+);7.1 数值微分【例7.27.2】生成一个5阶魔方矩阵，按列进行差分运算。M=magic(5)M=17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 9DM=diff(M)%计算V的一阶差分DM=6 19 6 6 1
4、 19 1 6 6 6 6 6 6 1 19 1 6 6 19 67.2 数值积分基本思想是将整个积分区间a，b分成n个子区间xi，xi+1，i=1，2，n，其中x1=a，xn+1=b。这样求定积分问题就分解为下面的求和问题：矩形法是用矩形面积近似曲边梯形的面积，如图7.2（a）所示；梯形法是用斜边梯形面积近似曲边梯形的面积，如图7.2（b）所示；而辛普生法是用抛物线近似曲边。7.2.1 数值积分的原理11()d()diinbxaxiSf xxf xx7.2 数值积分1梯形法2自适应辛普生法3高斯-克朗罗德（Gauss-Kronrod）法7.2.1 数值积分的原理7.2 数值积分1自适应积分算
5、法基于全局自适应积分算法的integral函数来求定积分。函数的调用格式为q=integral(fun,xmin,xmax)q=integral(fun,xmin,xmax,Name,Value)fun是被积函数，xmin和xmax分别是定积分的下限和上限。选项用于设置积分器的属性7.2.2 定积分的数值求解实现7.2 数值积分【例7.37.3】求（1）建立被积函数文件fe.m。function f=fe(x)f=x.*sin(x)./(1+abs(cos(x);（2）调用数值积分函数integral求定积分。q=integral(fe,0,pi)q=2.1776例7.3也可使用以下命令求解：
6、q=integral(x)(x.*sin(x)./(1+abs(cos(x),0,pi)7.2 数值积分 a=8755;b=6810;format long funLength=(x)sqrt(a2.*sin(x).2+b2.*cos(x).2);L=4*integral(funLength,0,pi/2)L=4.908996526868900e+047.2 数值积分2 2高斯高斯-克朗罗德法克朗罗德法quadgk函数求振荡函数的定积分。该函数的调用格式为q=quadgk(fun,xmin,xmax)q,errbnd=quadgk(fun,xmin,xmax,Name,Value)errbnd
7、返回近似误差范围。积分上下限可以是Inf或Inf，如果积分上下限是复数，则quadgk在复平面上求积分。7.2.2 定积分的数值求解实现7.2 数值积分【例7.57.5】求定积分 format long;I=quadgk(x)exp(x).*log(x),0,1)I=1.31790216241408110e ln dxx x7.2 数值积分3梯形积分法函数trapz对的离散数据用梯形法求定积分，调用格式如下。（1）T=trapz(Y)用于求均匀间距的积分。通常，Y Y是向量，采用单位间距（即间距为1），计算Y Y的近似积分。若Y Y是矩阵，则输出参数T是一个行向量，T的每个元素分别存储Y Y的
8、每一列的积分结果。（2）T=trapz(X,Y)用于求非均匀间距的积分。X X、Y Y满足函数关系Y=f(X)，按X指定的数据点间距，对Y求积分。7.2.2 定积分的数值求解实现7.2 数值积分【例7.67.6】从地面发射一枚火箭，表7.2记录了080秒火箭的加速度。试求火箭在第80秒时的速度。t=0:10:80;a=30.00,31.63,33.44,35.47,37.75,40.33,43.29,46.69,50.67;v=trapz(t,a)v=3.0893e+037.2 数值积分4 4累计梯形积分累计梯形积分在MATLAB中，提供了对数据积分逐步累计的函数cumtrapz。该函数调用格
9、式如下。（1）Z=cumtrapz(Y)用于求均匀间距的累计积分。（2）Z=cumtrapz(X,Y)用于求非均匀间距的累计积分。7.2.2 定积分的数值求解实现7.2 数值积分【例7.67.6】从地面发射一枚火箭，表7.2记录了080秒火箭的加速度。试求火箭在第80秒时的速度。计算例7.6中各时间点的速度，命令如下：v=cumtrapz(t,a)v=cumtrapz(t,a)v=1.0e+03*1 至 7 列 0 0.3081 0.6335 0.9780 1.3442 1.7346 2.1526 8 至 9 列 2.6026 3.08947.2 数值积分integral2integral2、
10、quad2dquad2d函数用于求二重函数用于求二重积分的积分的数值解，数值解，integral3integral3函数函数用于求三重积用于求三重积分的分的数值解。函数的调用格式为数值解。函数的调用格式为：q q=integral2(fun,xmin,xmax,ymin,ymax,Name,Value=integral2(fun,xmin,xmax,ymin,ymax,Name,Value)q q=quad2d(fun,xmin,xmax,ymin,ymax=quad2d(fun,xmin,xmax,ymin,ymax,),)q,errbnd=quad2d(fun,xmin,xmax,ymin
11、,ymax,Name,Valueq,errbnd=quad2d(fun,xmin,xmax,ymin,ymax,Name,Value)q=q=integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax)integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax)q q=integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax,Name,Value=integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax,Name,Value)fun为被积函数，xmin,xmax为x的积分区域，y
12、min,ymax为y的积分区域，zmin,zmax为z的积分区域，选项Name的用法及可取值与函数integral相同。输出参数q返回积分结果，errbnd用于返回计算误差。7.2.3 多重定积分的数值求解实现7.2 数值积分【例7.77.7】计算二重定积分 fxy=(x,y)exp(-x.2/2).*sin(x.2+y);I=integral2(fxy,-2,2,-1,1)I=1.5745212/2212esin()d dxxyx y【例7.87.8】计算三重定积分 fxyz=(x,y,z)4*x.*z.*exp(-z.*z.*y-x.*x);I3=integral3(fxyz,0,pi,0
13、,pi,0,1)I3=1.73282210004ed d dz y xxzx y z 7.3 离散傅立叶变换fft函数用于求离散傅立叶变换，调用格式如下。fft(X,n,dim)通常，X X是向量，若X X是矩阵，fft(X)应用于矩阵的每一列。选项n指定变换的数据点数，若X X的长度小于n，则不足部分补上零；若大于n，则删去超出n的那些元素。选项dim指定操作方向，dim=1时，该函数作用于X X的每一列；dim=2时，则作用于X X的每一行。7.2.2 定积分的数值求解实现7.3 离散傅立叶变换【例7.97.9】给定数学函数x(t)=3cos(2 20t+/3)+5sin(2 50t)，在
14、01s时间范围内采样128点，用fft作快速傅立叶变换，绘制相应的振幅频率图。N=128;%采样点数T=1;%采样时间终点t=linspace(0,T,N);x=3*cos(2*pi*20*t+pi/3)+5*sin(2*pi*50*t);%求各采样点样本值dt=t(2)-t(1);%采样周期fs=1/dt;%采样频率(Hz)X=fft(x);%计算x的快速傅立叶变换F=X(1:N);f=fs/N*(0:N-1);%得到信号的频率成分plot(f,abs(F),-*)%绘制振幅频率图xlabel(Frequency);ylabel(|F(k)|)7.4 常微分方程的数值求解4阶龙格-库塔（Ru
15、nge-Kutta）公式为7.4.1 龙格-库塔法简介12112312241311234(,)(,)2(,)2(,)(22)6iiiiiiiiiikf x yhkf xykhkf xykkf xyhkhyykkkk7.4 常微分方程的数值求解1.求解函数求解常微分方程数值解的函数用法相同，其基本调用格式为T,Y=solver(odefun,tspan,y0,options)T,Y,TE,YE,IE=solver(odefun,tspan,y0,options)其中，solver是根据待求解问题的性质选用的求解函数。7.4.2 常微分方程数值求解的实现7.4 常微分方程的数值求解2.求解函数的参
16、数（1）输出参数T为列向量，返回时间点；Y为矩阵，Y的每一行向量返回T中各时间点相应的状态。TE为事件发生的时刻，YE为对应时刻的解，IE为触发事件的索引。（2）输入参数odefun为函数文件名或匿名函数，函数文件头通常采用f(t，y)的形式，即形参t为时间参量，形参y为待求解问题的自变量。tspan指定求解区间，用二元向量t0 tf 表示。y0是初始状态列向量。options用于设置积分求解过程和结果的属性。7.4.2 常微分方程数值求解的实现7.4 常微分方程的数值求解【例7.107.10】设有初值问题试求其数值解，并与精确解相比较，精确解为（1）建立函数文件funst.m。functio
17、n yp=funst(x,y)yp=y.*tan(x)+sec(x);（2）求解微分方程。y0=pi/2;x,y=ode23(funst,0,1,y0);%求数值解 y1=(x+pi/2)./cos(x);%求精确解 x,y,y1ans=0 1.5708 1.5708 0.1000 1.6792 1.6792 0.2000 1.8068 1.8068 0.3000 1.9583 1.9583 0.4000 2.1397 2.1397 0.5000 2.3596 2.3597 0.6000 2.6301 2.6302 0.7000 2.9689 2.9690 0.8000 3.4027 3.40
18、29 0.9000 3.9745 3.9748 1.0000 4.7573 4.75810tgsec,01,2xyyxxxy 7.4 常微分方程的数值求解【例7.117.11】求描述振荡器的Van der Pol方程：令x1=y，x2=y，则可写出Van der Pol方程的状态方程：（1）建立函数文件verderpol.m。function xprime=verderpol(t,x)global mu;xprime=x(2);mu*(1-x(1)2)*x(2)-x(1);（2）求解微分方程。global mu;mu=2;y0=1;0;t,x=ode45(verderpol,0,20,y0);
19、subplot(1,2,1);plot(t,x);%绘制系统时间响应曲线 subplot(1,2,2);plot(x(:,1),x(:,2)%绘制系统相平面曲线7.4 常微分方程的数值求解【例7.127.12】某非刚性物体的运动方程为其初始条件为x(0)=0，y(0)=0，z(0)=。取=8/3，=28，=10，试绘制系统相平面图。将运动方程写成矩阵形式，得（1）建立模型的函数文件lorenz.m。function xdot=lorenz(t,x)xdot=-8/3,0,x(2);0,-10,10;-x(2),28,-1*x;（2）解微分方程组。t,x=ode23(lorenz,0,80,0;0;eps);plot3(x(:,1),x(:,2),x(:,3);axis(0,50,-20,30,-35,28);
