第六章 matlab数值计算.ppt

《第六章 matlab数值计算.ppt》由会员分享，可在线阅读，更多相关《第六章 matlab数值计算.ppt（56页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、Matlab 仿真及其应用主讲：陈孝敬e-mail:主要内容：主要内容：多项式运算；多项式运算；插值运算；插值运算；数据分析；数据分析；微分方程组数值解；微分方程组数值解；第第6 6章章 数值计算数值计算6.16.1多项式运算多项式运算1 1多项式表示法多项式表示法matlabmatlab语言把多项式表达成一个行向量，语言把多项式表达成一个行向量，该向量中的元素是按多项式降幂排列的。该向量中的元素是按多项式降幂排列的。f(x)=af(x)=an nx xn n+a+an-1n-1x xn-1n-1+a+a0 0 可用行向量可用行向量 p=ap=an n a an-1n-1 a a1 1+a+a

2、0 0 表示表示 poly poly 产生特征多项式系数向量产生特征多项式系数向量 特征多项式一定是特征多项式一定是n+1n+1维的维的 特征多项式第一个元素一定是特征多项式第一个元素一定是1 1例例6-16-1:a a=1 2 3;4 5 6;7 8 0;=1 2 3;4 5 6;7 8 0;p=poly(a)p=poly(a)p=1.00 -6.00 -72.00 -27.00p=1.00 -6.00 -72.00 -27.00p p是多项式是多项式p(x)=xp(x)=x3 3-6x-6x2 2-72x-27-72x-27的的MatlabMatlab描述方法，我们描述方法，我们可用：可用

3、：p1=poly2str(p,p1=poly2str(p,x x)函数文件，显示函数文件，显示数学多项式的形式数学多项式的形式p1=x3-6 x2-72 x p1=x3-6 x2-72 x 27 272 2多项式求值多项式求值MatlabMatlab用用polyval()polyval()来计算多项式的值来计算多项式的值y=polyval(p,x),py=polyval(p,x),p为多项式系数行向量，为多项式系数行向量，x x代入多项式的值。代入多项式的值。Y=polyvalm(p,X),Y=polyvalm(p,X),把矩阵把矩阵X X代入多项式代入多项式p p中进行计算，其中矩阵中进行计

4、算，其中矩阵X X必须是方阵。必须是方阵。例例6-26-2:p:p=2=2 2 2 3;3;A=1 0;0 2;A=1 0;0 2;r_a=polyval(p,A);r_a=polyval(p,A);r_b=polyvalm(p,A)r_b=polyvalm(p,A)2 2多项式乘法和除法多项式乘法和除法在在MatlabMatlab中用函数中用函数conv()conv()和和deconv()deconv()进行多项式乘法和除法。进行多项式乘法和除法。例例6-36-3:p1p1=1=1 2 2 3;3;p2=6 0 0 3 3 4;p2=6 0 0 3 3 4;w=conv(p1,p2);w=c

5、onv(p1,p2);q,r=deconv(p2,p1);q,r=deconv(p2,p1);sq=poly2str(q,sq=poly2str(q,x x););sr=poly2str(r,sr=poly2str(r,x x););3 3多项式的微积分多项式的微积分在在MatlabMatlab中用函数中用函数polyder()polyder()来计算多项式的微分。来计算多项式的微分。用函数用函数polyint()polyint()返回多项式的不定积分返回多项式的不定积分例例6-46-4:p1:p1=1=1 2 2 3;3;p2=6 0 0 3 3 4;p2=6 0 0 3 3 4;k1=po

6、lyder(p1);k1=polyder(p1);k2=polyder(p1,p2);k2=polyder(p1,p2);k3,k4=polyder(p2,p1);k3,k4=polyder(p2,p1);sq=poly2str(k3,sq=poly2str(k3,x x););sr=poly2str(k4,sr=poly2str(k4,x x););例例6-56-5:p p=1=1 2 2 3;3;k1=ployint(p1);k1=ployint(p1);k2=ployint(p1,2);k2=ployint(p1,2);3 3多项式的微积分多项式的微积分在在MatlabMatlab中用函

7、数中用函数polyder()polyder()来计算多项式的微分。来计算多项式的微分。用函数用函数polyint()polyint()返回多项式的不定积分返回多项式的不定积分4.roots 求多项式的根求多项式的根例例6-66-6:a=1 2 3;4 5 6;7 8 0;p=poly(a)p=1.00 -6.00 -72.00 -27.00r=roots(p)r=12.12 -5.73 显然显然 r是矩阵是矩阵a的特征值的特征值 -0.39 例例6-76-7:用函数用函数roots（）来计算多项式（）来计算多项式x5+x4-8x3-6x2+8x+24的根的根p=1 1-8-6 8 24;sp=

8、poly2str(p,x);r=roots(p)例例6-8 由根由根-2 2 1创建多项式。创建多项式。r=-2 2 1;p=poly(r);sp=poly2str(p,x);由根创建多项式由根创建多项式,在在Matlab中用函数中用函数poly()实现实现5.多项式曲线拟合多项式曲线拟合在在Matlab中用函数中用函数polyfit实现采用最小二乘法对给定数据实现采用最小二乘法对给定数据进行多显示拟合，进行多显示拟合，p=ployfit(x,y,n)，采用，采用n次多项式次多项式p来拟来拟合数据线合数据线x和和y。例例6-9 用函数用函数polyfit()对给定数据进行多项式拟合。对给定数据

9、进行多项式拟合。x=0:0.2:10;y=0.25*x+20*sin(x);%5阶多项式拟合阶多项式拟合p1=polyfit(x,y,5);y1=polyval(p1,x);%8阶多项式拟合阶多项式拟合p2=polyfit(x,y,8);y2=polyval(p2,x);hold on plot(x,y,ro);plot(x,y1,r-);plot(x,y2,b:);xlabel(x);ylabel(y);legend(原始数据原始数据,5阶多项式拟合阶多项式拟合,8阶多项式拟合阶多项式拟合);6.26.2 插值运算插值运算插值是根据已知输入插值是根据已知输入/输出数据集和当前输入估计输出值，

10、输出数据集和当前输入估计输出值，其中，当前输入应在已知数据集内。在信号和图像处理其中，当前输入应在已知数据集内。在信号和图像处理中是经常用到的。中是经常用到的。Matlab提供了一维、二维、提供了一维、二维、三次样条等许多插值选择三次样条等许多插值选择插值方法Matlab函数插值方法Matlab函数一维插值interp1使用FFT方法的一维插值interpft快速一维插值interpq分段三次Hermite插值pchip二维插值interp2三次样条插值spline三维插值interp3N维插值interpn6.2.16.2.1 一维插值一维插值一维插值原理请参照一维插值原理请参照P115，在

11、，在Matlab中用中用interp1()来实来实现一维插值。现一维插值。调用格式：调用格式：yiyiinterp1(x,y,xi)interp1(x,y,xi)已知数据向量（已知数据向量（x,yx,y），计算并返），计算并返回在插值向量回在插值向量xixi处的函数值处的函数值yi=interp1(x,y,xi,yi=interp1(x,y,xi,methodmethod)yi=interp1(x,y,xi,yi=interp1(x,y,xi,methodmethod,extrapextrap)methodmethod用于指定插值算法，其值可以是：用于指定插值算法，其值可以是：nearestn

12、earest最近插值最近插值linearlinear线性插值（默认值）线性插值（默认值）splinespline分段三次样条插值分段三次样条插值pchippchip分段三次分段三次HermiteHermite插值插值cubiccubic与与pchippchip相同相同一维插值一维插值方法比较方法比较例例6-10 用不同插值法对数据进行一维插值用不同插值法对数据进行一维插值x=0:1:10;x=0:1:10;y=0 0.8415 0.9093 0.1411 -0.7568-0.9589-0.2794 0.657 y=0 0.8415 0.9093 0.1411 -0.7568-0.9589-0.

13、2794 0.657 0.9894 0.4121 -0.5440;0.9894 0.4121 -0.5440;plot(x,y,co),hold onplot(x,y,co),hold onfplot(sin,0 10)fplot(sin,0 10)xi=0:0.15:10;xi=0:0.15:10;yi=interp1(x,y,xi);yi=interp1(x,y,xi);figure,plot(xi,yi,r+),text(0.7028,0.4649,figure,plot(xi,yi,r+),text(0.7028,0.4649,线性插值线性插值rightarrow)rightarrow

14、)yi2=interp1(x,y,xi,nearst);yi2=interp1(x,y,xi,nearst);figure,plot(xi,yi2,c*),text(3.537,0.1374,leftarrowfigure,plot(xi,yi2,c*),text(3.537,0.1374,leftarrow最近插值最近插值)yi3=interp1(x,y,xi,cubic);yi3=interp1(x,y,xi,cubic);plot(xi,yi3,md),text(2.408,0.8333,leftarrowplot(xi,yi3,md),text(2.408,0.8333,leftarr

15、ow三次插值三次插值)yi4=interp1(x,y,xi,spline);yi4=interp1(x,y,xi,spline);figure,plot(xi,yi4,kh),text(4.62,0.8158,figure,plot(xi,yi4,kh),text(4.62,0.8158,三次样条插值三次样条插值rightarrow)rightarrow)初始数据对于插值的影响初始数据对于插值的影响x=0:2:10;y=sin(x);plot(x,y,go),hold onezplot(sin,0 10)xi=0:0.15:10;yi=interp1(x,y,xi);figure,plot(x

16、i,yi,r+),text(0.5876,0.2537,leftarrow线性插值线性插值)yi2=interp1(x,y,xi,nearst);figure,plot(xi,yi2,c*),text(6.947,-0.258,leftarrow最近插值最近插值)yi3=interp1(x,y,xi,pchip);figure,plot(xi,yi3,md),text(2.408,0.8333,leftarrow三次插值三次插值)yi4=interp1(x,y,xi,spline);figure,plot(xi,yi4,kh),text(1.601,1.138,leftarrow三次样条插值三

17、次样条插值)函数函数spline与与pchipspline()的调用格式为：的调用格式为：yi=spline(x,y,xi)此函数等同于此函数等同于yi=interp1(x,y,xi,spline)pp=spline(x,y)返回三次样条插值的分段多项式形式的向量返回三次样条插值的分段多项式形式的向量spline函数可以保证插值函数的三阶导数连续函数可以保证插值函数的三阶导数连续pchip()的调用格式为：的调用格式为：yi=pchip(x,y,xi)此函数等同于此函数等同于yi=interp1(x,y,xi,pchip)pp=pchip(x,y)返回三次样条插值的分段多项式形式的向量返回三次

18、样条插值的分段多项式形式的向量函数函数pchip与与spline的区别的区别pchip与与spline的区别：三次样条在相邻的节点上并不保证单的区别：三次样条在相邻的节点上并不保证单调性；而调性；而Hermite分段三次样条则可保证插值的局部单调性分段三次样条则可保证插值的局部单调性例例6-11 利用点利用点(x=sin(k/6),y=cos(k/6)，其中，其中k=0 1 2 3来逼来逼近单位圆的前四分之一圆周。比较近单位圆的前四分之一圆周。比较pchip与与spline的差别。的差别。t=linspace(0,pi/2,4)x=cos(t);y=sin(t);xx=linspace(0,1

19、,40);plot(x,y,s,xx,pchip(x,y,xx);spline(x,y,xx)grid on,axis equallegend(Orignal data,spline,pchip)6.2.26.2.2 二维插值二维插值Matlab 中用中用interp2函数实现函数实现调用格式：调用格式：zi=interp2(x,y,z,xi,yi,method)method算法属性值可以是；算法属性值可以是；nearest最近插值最近插值linear线性插值（默认）线性插值（默认）spline三次样条插值三次样条插值(spline)cubic立方插值立方插值例例6-12 假设有一组分度系数的

20、假设有一组分度系数的“海底深度测量数据海底深度测量数据”，由以，由以下一段程序生成：下一段程序生成：randn(state,2);x=-5:5;y=-5:5;X,Y=meshgrid(x,y);Z=-500+1.2*exp(-(X-1).2+(Y-2).2)-0.7*exp(-(exp(X+2).2+(Y+1).2);surf(X,Y,Z)，view(-25,25)试由插值方式绘制海底形状图。试由插值方式绘制海底形状图。xi=linspace(-5,5,50);yi=linspace(-5,5,50);XI,YI=meshgrid(xi,yi);ZI=interp2(X,Y,Z,XI,YI,*

21、cubic);surf(XI,YI,ZI)view(-25,25)例例6-13 不同插值方法对数据进行二维插值：不同插值方法对数据进行二维插值：x,y=meshgrid(-3:0.8:3);z=peaks(x,y);xi,yi=meshgrid(-1:0.25:3);zi_nearest=interp2(x,y,z,xi,yi,nearset);zi_linear=interp2(x,y,z,xi,yi);zi_spline=interp2(x,y,z,xi,yi,spline);figure;hold on;subplot(2,2,1);meshc(x,y,z);title(原始数据原始数据

22、);subplot(2,2,2);meshc(xi,yi,zi_nearest);title(最邻近插值最邻近插值);subplot(2,2,3);meshc(xi,yi,zi_linear);title(线性插值线性插值);subplot(2,2,4);meshc(xi,yi,zi_spline);title(三次样条插值三次样条插值);6.36.3 数据分析数据分析6.3.1 基本数据分析函数请参照基本数据分析函数请参照P118-119表表6-3例例6-14 应用求最大值、最小值、平均值、中间值、元素和等应用求最大值、最小值、平均值、中间值、元素和等函数。函数。x=1:40;y=randn

23、(1,40);figure;hold on;plot(x,y);y_max,I_max=max(y);plot(x(I_max),y_max,o);y_min,I_min=min(y);plot(x(I_min),y_min,*);xlabel(x);xlabel(y);legend(原始数据原始数据,最大值最大值,最小值最小值);y_mean=mean(y);y_median=median(y);y_sum=sum(y);例例6-15 计算向量的标准差和方差计算向量的标准差和方差x=1:10;mean_x=mean(x);r=0;for i=1:10 r=r+(x(i)-mean_x)2;e

24、ndr1=sqrt(r/10)r2=sqrt(r/9)r3=std(x);r4=r1*2r5=r22r6=var(x)r7=var(x,1)例例6-16 对向量排序对向量排序a=4 7-5;b=1 0-5 10-6;p=roots(a);b,c=sort(a)b,c=sort(p)例例6-17 对字符串排序对字符串排序a=hello;world;hally;clayt;Daney;b=sortrows(a)%字符串首先按第一个字母排序；若第一个字字符串首先按第一个字母排序；若第一个字%符相同，再按第二个字符排序，以此类推符相同，再按第二个字符排序，以此类推c=sortrows(a,2)%字符串

25、首先按第二个字母排序；若第二字符串首先按第二个字母排序；若第二%个字个字%符相同，再按第三个字符排序，以此类推符相同，再按第三个字符排序，以此类推例例6-18 计算随机变量样本的协方差和相关系数矩阵计算随机变量样本的协方差和相关系数矩阵a=rand(10,4);c1=cov(x);c2=cov(x,1);c3=corrcoef(x);6.3.2 协方差和相关系数矩阵协方差和相关系数矩阵若给定若给定n个个m维随机变量样维随机变量样 ，定义如下矩，定义如下矩阵为其协方差矩阵：阵为其协方差矩阵：其中，其中，同时定义，同时定义矩阵为相关系数矩阵矩阵为相关系数矩阵 其中其中6.3.3 有限差分和梯度有限

26、差分和梯度在在Matlab中用函数中用函数diff()来计算差分，用函数来计算差分，用函数gradient()来计来计算梯度算梯度例例6-19 利用有限差分近似计算正弦函数的导数利用有限差分近似计算正弦函数的导数x=1:0.1:10;y=sin(x);y_der=diff(y)./diff(x);hold on;x_der=x(1:(end-1);plot(x,y,b-);plot(x_der,y_der,b-.);axis(0 10-1.2 1.4);legend(正弦函数正弦函数,正弦函数导数正弦函数导数);对于向量对于向量X，diff(X)表示了表示了X(2)-X(1)X(3)-X(2)

27、.X(n)-X(n-1).对于矩阵对于矩阵X，diff(X)表示了表示了X(2:n,:)-X(1:n-1,:)diff(x,n,dim)得到矩阵得到矩阵x在在dim维上的维上的n阶差值阶差值例例6-20二维高斯函数的梯度场二维高斯函数的梯度场v=-2:0.25:2;x,y=meshgrid(v,v);%产生字变量产生字变量x，yz=exp(-(x.2+y.2+0.5*x.*y);%二维高斯函数二维高斯函数px py=gradient(z,0.25);%梯度场梯度场quiver(v,v,px,py);%画出梯度场画出梯度场6.3.4 信号滤波和卷积信号滤波和卷积函数名功能描述插值方法Matlab

28、函数filter一维数字滤波器convNN维卷积filter2二维数字滤波器deconv反卷积和多项式除法conv一维卷积和多项式乘法detrend去除信号中的直流或者线性conv2二维卷积例例6-216-21 对一带噪音的正弦信号进行对一带噪音的正弦信号进行5 5阶平均值滤波阶平均值滤波t=0:0.1:10;%时间时间 n=6*randn(size(t);%高斯白噪音高斯白噪音s=40*sin(t);%无噪音的正弦信号无噪音的正弦信号x=40*sin(t)+n;%在正弦信号中添加噪音在正弦信号中添加噪音a=1;b=1/5 1/5 1/5 1/5 1/5;%滤波器滤波器y=filter(b,a

29、,x);%滤波滤波plot(t,s,g-.);%画无噪音信号画无噪音信号hold on;plot(t,x,b-);%画染噪信号画染噪信号 plot(t,y,r:);%画滤波后的信号画滤波后的信号 axis(0 10-65 65);xlabel(时间时间(s);legend(无噪信号无噪信号,有噪信号有噪信号,滤波信号滤波信号);例例6-226-22 计算向量的卷积计算向量的卷积u=ones(1,15);%阶跃信号阶跃信号v=1/20:1/20:1;%线性信号线性信号w=conv(u,v);%卷积卷积figure;subplot(3,1,1);stem(u);title(u);subplot(3

30、,1,2);stem(v);title(v);subplot(3,1,3);stem(w);title(w);6.3.5 傅立叶变换傅立叶变换离散傅立叶变换的实现离散傅立叶变换的实现离散傅立叶变换函数如表离散傅立叶变换函数如表:函数名功能描述插值方法Matlab函数fft一维离散快速傅立叶变换ifft一维离散快速傅立叶逆变换fft2二维离散快速傅立叶变换ifft2二维离散快速傅立叶逆变换fftnn维离散快速傅立叶变换ifftnn维离散快速傅立叶逆变换例例6-23 给定数学函数给定数学函数x(t)=12sin(210t+/4)+5cos(240t)取取N=128，试对，试对t从从01秒采样，用秒

31、采样，用fft作快速傅立叶变换，绘作快速傅立叶变换，绘制相应的振幅制相应的振幅-频率图。频率图。在在01秒时间范围内采样秒时间范围内采样128点，从而可以确定采样周期和采点，从而可以确定采样周期和采样频率。由于离散傅立叶变换时的下标应是从样频率。由于离散傅立叶变换时的下标应是从0到到N-1，故，故在实际应用时下标应该前移在实际应用时下标应该前移1。又考虑到对离散傅立叶变。又考虑到对离散傅立叶变换来说，其振幅换来说，其振幅|F(k)|是关于是关于N/2对称的，故只须使对称的，故只须使k从从0到到N/2即可。即可。程序如下：程序如下：N=128;%采样点数采样点数T=1;%采样时间终点采样时间终点

32、t=linspace(0,T,N);%给出给出N个采样时间个采样时间ti(I=1:N)x=12*sin(2*pi*10*t+pi/4)+5*cos(2*pi*40*t);%求各采样点求各采样点样本值样本值xdt=t(2)-t(1);%采样周期采样周期f=1/dt;%采样频率采样频率(Hz)X=fft(x);%计算计算x的快速傅立叶变换的快速傅立叶变换XF=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),-*)%绘制振幅绘制振幅-频率图频率图xlabel(Frequency);ylabel

33、(|F(k)|)例例6-24分析一图像的频谱分析一图像的频谱load imdemos saturn2imshow(saturn2)B=fftshift(fft2(saturn2);imshow(log(abs(B),notruesize)6.46.4 功能函数功能函数功能函数就是可将其他函数作为输入变量的函数功能函数就是可将其他函数作为输入变量的函数1.函数的表示函数的表示在在Matlab中，函数可以通过中，函数可以通过M文件、匿函数和函数文件、匿函数和函数inline()来表示来表示例例6-25 用上述用上述3种方式表示函数种方式表示函数y=sin(x)+e-2xfunction y=fun

34、express(x)y=sin(x)+exp(-2*x);fh=(x)sin(x)+exp(-2*x);g=inline(sin(x)+exp(-2*x);2.函数画图函数画图例例6-26 绘制函数绘制函数 在区间在区间0，10上的图像上的图像 f=(x)cos(x+1)./(x.2+1);fplot(f,-5 5,1e-4,r-);title(函数函数y=cos(x+1)/(x2+1);xlabel(x);ylabel(y);grid;3.函数最小值和零点函数最小值和零点求解一元函数的最小值求解一元函数的最小值可以通过函数可以通过函数fminbndfminbnd来求一元函数来求一元函数y y

35、=f f(x x)在指定区间在指定区间 a a,b b 上上的函数局部极小值的函数局部极小值,该函数返回函数在极小值点时自变量该函数返回函数在极小值点时自变量x x的值的值,调用格式为调用格式为:x=fminbnd(:x=fminbnd(fun,a,b).fun,a,b).例例6-276-27.求求humpshumps函数在开区间函数在开区间(0.3,1)(0.3,1)内的最小值内的最小值.humps.humps是是MatlabMatlab内置的内置的M M文件函数文件函数,实际上是实际上是y=1/(x-0.3)2+0.01)+1/(x-0.9)2+0.04)-6.y=1/(x-0.3)2+0

36、.01)+1/(x-0.9)2+0.04)-6.x=fminbnd(humps,0.3,1)求解多元函数的最小值求解多元函数的最小值函数函数fminsearchfminsearch用于求多元函数在向量用于求多元函数在向量x0 x0附近的最小值附近的最小值.它指定它指定一个开始的向量一个开始的向量(x0),(x0),并非指定一个区间并非指定一个区间.此函数返回一个向量此函数返回一个向量,为此多元函数局部最小函数值对应的自变量的取值为此多元函数局部最小函数值对应的自变量的取值,调用格式为调用格式为x=fminsearch(fun,x0)x=fminsearch(fun,x0)例例6-286-28.

37、把一个把一个3 3个自变量的函数创建在一个个自变量的函数创建在一个M M文件里，文件里，求这个函数在求这个函数在1,-1,01,-1,0点附近的最小值可以得到点附近的最小值可以得到:.function b=three(v)x=v(1);y=v(2);z=v(3);b=x*x+2.5*sin(y)-z*z*x*y*y;v=1-1 0;fminsearch(three,v)求一元函数的零点求一元函数的零点函数函数fzerofzero用于求一元函数的零点，用于求一元函数的零点，x=zero(x=zero(funfun,x0),x0),在在x0 x0点附近寻找函数的零点点附近寻找函数的零点x=zero

38、(x=zero(funfun,x0,x0，x1),x1),在在x0,x1x0,x1区间内寻找寻找函数的区间内寻找寻找函数的零点零点x=zero(x=zero(funfun,x0,options),x0,options),用用optionsoptions参数指定寻找零点的参数指定寻找零点的优化器参数。优化器参数。定义：若函数定义：若函数y yf f（x x）的定义域为）的定义域为D D，x0Dx0D，且，且f f（x0 x0）0 0，则称则称x xx0 x0是函数是函数y yf f（x x）的零点。）的零点。零点的几何意义，函数图象与零点的几何意义，函数图象与x x轴的公共点（交点或切点）的横轴

39、的公共点（交点或切点）的横坐标。坐标。从方程看，零点就是方程从方程看，零点就是方程f f（x x）0 0在在D D上的根上的根例例6-296-29 求函数求函数 的零点的零点f=(x)1./(x+4).2+1)+1./(x-4).2+1)-0.5;fplot(f,-10 10);xlabel(x);ylabel(f(x);x1=fzero(f,-5)x2=fzero(f,-3)x3=fzero(f,3)x4=fzero(f,5)x1_region=fzero(f,-10,-3)x2_region=fzero(f,-3,0)x3_region=fzero(f,0,3)x4_region=fzer

40、o(f,3,10)例例6-306-30.显示显示fminbnd()fminbnd()的寻找最小值的过程的寻找最小值的过程.f=inline(sin(x)+3);%用内联函数表达用内联函数表达x=fminbnd(f,2,5)%如果需要显示计算过程，可以如果需要显示计算过程，可以optimset()来设定参数来设定参数X=fminbnd(f,2,5,optimset(display,iter);优化器参数优化器参数在求一元函数、多元函数最小值以及一元函数零点时，都可以在求一元函数、多元函数最小值以及一元函数零点时，都可以设定求优化器的参数，参数优化器设定求优化器的参数，参数优化器optimset(

41、)optimset()4.数值积分数值积分数值积分在数值计算中有着重要作用数值积分在数值计算中有着重要作用,许多数值计算问题许多数值计算问题可以转化为数值积分问题可以转化为数值积分问题,如常微分方程初值问题等如常微分方程初值问题等通常可用逼近多项式通常可用逼近多项式Pn(x)来代替被积函数来代替被积函数f(x)，计算积，计算积分分构造数值积分的方法很多，主要有构造数值积分的方法很多，主要有Newton-Cotes系列数系列数值积分法、值积分法、Gauss积分法和积分法和Romberg积分法等积分法等Matlab函数公式quad自适应Simpson求积公式（低阶）quadl自适应Lobatto求

42、积公式；精度高，最常用trapz梯形求积公式；速度快，精度差cumtrapz梯形法求一个区间上的积分曲线cumsum等宽距法求一个区间上的积分曲线，精度很差fnint利用样条函数求不定积分；与spline，ppval配合使用，主要应用于表格“函数”积分例例6-316-31.比较不同的方法，求比较不同的方法，求sinsin在区间在区间0 0 上的积分上的积分.y1=(x)sin(x);t=0:pi/100:pi;y2=sin(t);q1=quad(y1,0,pi);q2=quadl(y1,0,pi);q3=trapz(0:pi/100:pi,y2);q4=cumtrapz(0:pi/100:pi

43、,y2);例例6-326-32 求求 .fun=inline(1./(sqrt(x).*(exp(x)+1);quadl(fun,0,1)quadl(fun,eps,1)广义积分广义积分奇点积分奇点积分 无穷积分无穷积分 例例6-336-33 求求 .可先选取一个有限的积分区间，如可先选取一个有限的积分区间，如0,100计算；在选择一个较大计算；在选择一个较大的积分区间，如的积分区间，如0 200计算，如两次计算结果的差满足一定的精计算，如两次计算结果的差满足一定的精度要求，则可认为此值即为无穷积分的值度要求，则可认为此值即为无穷积分的值多重数值积分多重数值积分二重积分函数：二重积分函数：SS

44、=dblquad(fun,xmin,xmax,ymin,ymax,tol,method)三重积分函数：三重积分函数：SSS=triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol,method)说明说明:MatlabMatlab只能处理积分限为常数的多重积分，对内积分上下限为外只能处理积分限为常数的多重积分，对内积分上下限为外积分变量的积分问题积分变量的积分问题,需自行编写函数求解。需自行编写函数求解。例例6-346-34.求二维高斯函数在矩形区间求二维高斯函数在矩形区间-1 1-1 1-1 1-1 1和圆形区域和圆形区域.上的二重积分。上的二重积分。

45、f=(x,y)1/sqrt(pi)*exp(-x.2)*1/sqrt(pi)*exp(-y.2);dblquad(f,-1,1,-1,1,1e-6,quadl)6.56.5 微分方程组数值解微分方程组数值解Matlab常微分方程求解问题分类（详细见P139）初值问题：定解附加条件在自变量的一端 一般形式为：初值问题的数值解法一般采用步进法，如Runge-Kutta法边值问题：在自变量两端均给定附加条件一般形式：边值问题可能有解、也可能无解，可能有唯一解、也可能有无数解边值问题有3种基本解法迭加法打靶法松弛法1.Matlab求解常微分方程初值问题方法将待求解转化为标准形式，并“翻译”成Matla

46、b可以理解的语言，即编写odefile文件选择合适的解算指令求解问题根据求解问题的要求，设置解算指令的调用格式Matlab求解初值问题函数指令含义指令含义解算ode23普通2-3阶法解ODEodefileODE文件格式ode45普通4-5阶法解ODE选项odeset创建、更改ODE选项的设置ode113普通变阶法解ODEodeget读取ODE选项的设置ode23t解适度刚性ODE输出odeplotODE的输出时间序列图ode15s变阶法解刚性ODEodephas2ODE的二维相平面图ode23s低阶法解刚性ODEodephas3ODE的三维相空间图ode23tb低阶法解刚性ODEodeprin

47、t在Matlab指令窗显示结果例例6-356-35 求求 的解的解 。首先，将上述方程改写成如下一阶常微分方程组：首先，将上述方程改写成如下一阶常微分方程组：function dydt=ivpodefun(t,y)dydt=zeros(2,1);dydt(1)=y(2);dydt(2)=(1-y(1)2)*y(2)-y(1);t,y=ode45(ivpodefun,0 20,2;0);plot(t,y(:,1),-,t,y(:,2),-);title(常微分方程的解);xlabel(t);ylabel(y);legend(y,y的一阶导数);2.常微分方程的边界问题常微分方程的边界问题Matl

48、ab求解边值问题方法：求解边值问题方法：bvp4c函数函数把待解的问题转化为标准边值问题把待解的问题转化为标准边值问题 因为边值问题可以多解，所以需要为期望解指定一个初始猜测解。该因为边值问题可以多解，所以需要为期望解指定一个初始猜测解。该猜测解网（猜测解网（Mesh）包括区间）包括区间a,b内的一组网点（内的一组网点（Mesh points）和网）和网点上的解点上的解S(x)根据原微分方程构造残差函数根据原微分方程构造残差函数利用原微分方程和边界条件，借助迭代不断产生新的利用原微分方程和边界条件，借助迭代不断产生新的S(x)，使残差不，使残差不断减小，从而获得满足精度要求的解断减小，从而获得

49、满足精度要求的解 MatlabMatlab求解边值问题求解边值问题的基本指令的基本指令solinitsolinitbvpinitbvpinit（x,v,parametersx,v,parameters）生成生成bvp4cbvp4c调用指令所必须的调用指令所必须的“解猜测网解猜测网”solsolbvp4cbvp4c（odefunodefun，bcfunbcfun，solinitsolinit，optionsoptions，p1p1，p2p2，）给出微分方程边值问题的近似解给出微分方程边值问题的近似解 sxintsxintdevaldeval（solsol，xintxint）计算微分方程积分区间内

50、任何一点的解值计算微分方程积分区间内任何一点的解值初始解生成函数：初始解生成函数：bvpinit()bvpinit()solinitbvpinit（x,v,parameters）x指定边界区间指定边界区间a,b上的初始网络，通常是等距排列的上的初始网络，通常是等距排列的（1M）一维数组。注意：使）一维数组。注意：使x(1)=a，x(end)=b；格点要单调；格点要单调排列。排列。v是对解的初始猜测是对解的初始猜测solinit（可以取别的任意名）是（可以取别的任意名）是“解猜测网（解猜测网（Mesh）”。它是。它是一个结构体，带如下两个域：一个结构体，带如下两个域：solinit.x是表示初始