第5章-数值计算课件

第五章第五章 数值计算数值计算目 录 5.1 5.1 多项式运算多项式运算5.2 5.2 插插 值值 运运 算算5.3 5.3 数数 据据 分分 析析5.4 5.4 功功 能能 函函 数数 5.5 5.5 微分方程组数值解微分方程组数值解 习习 题题5.1 5.1 多项式运算多项式运算1.1.多项式多项式MatlabMatlab表示法表示法2.2.多项式求值多项式求值3.3.多项式乘法和除法多项式乘法和除法 4.4.多项式的微积分多项式的微积分5.5.多项式的根和由根创建多项式多项式的根和由根创建多项式6.6.多项式部分分式展开多项式部分分式展开 7.7.多项式曲线拟合多项式曲线拟合8.8.多项式曲线拟合图形用户接口多项式曲线拟合图形用户接口MatlabMatlab提供了下列关于多项式的函数：提供了下列关于多项式的函数：1.1.多项式表示法多项式表示法lMATLABMATLAB采采用用行行向向量量表表示示多多项项式式系系数数，多多项项式系数按降幂排列。式系数按降幂排列。lpoly2strpoly2str函函数数将将多多项项式式系系数数向向量量转转换换为为完完整形式。整形式。poly2str(1 0 2,s)ans=s2+22.2.多项式求值多项式求值polyvalpolyval函数计算多项式的值，其用法为：函数计算多项式的值，其用法为：l Y=polyval(P,X),P多项式系数行向量X代入多项式的值其意义为：当当P P为长度为为长度为n+1n+1的行向量的行向量 pp1 1,p,p2,2,.,p.,pn n,p,pn+1n+1 时，时，Y=Y=p p1 1x xn n+p p2 2x xn-1n-1+.+.+p pn nx+x+p pn+1n+1例子：P=1:6,poly2str(P,x),x=5,polyval(P,x)P=1 2 3 4 5 6ans=x5+2 x4+3 x3+4 x2+5 x+6x=5ans=4881例子：P=1:6,poly2str(P,x),x=2 3;4 5,polyval(P,x)P=1 2 3 4 5 6ans=x5+2 x4+3 x3+4 x2+5 x+6x=2 3 4 5ans=120 543 1818 4881lY=polyval(P,X),把矩阵或向量X中的每个元素逐个代入多项式中进行计算；lY Y=polyvalm(P,X)polyvalm(P,X)，把把矩矩阵阵X X作作为为整整体体代代入入多多项项式式中进行计算，中进行计算，X X必须为方阵。必须为方阵。例子：P=1:6,poly2str(P,x),x=2 3;4 5,polyval(P,x)P=1 2 3 4 5 6ans=x5+2 x4+3 x3+4 x2+5 x+6x=2 3 4 5ans=120 543 1818 4881例子：P=1:6,poly2str(P,x),x=2 3;4 5,polyvalm(P,X)P=1 2 3 4 5 6ans=x5+2 x4+3 x3+4 x2+5 x+6x=2 3 4 5ans=8256 10881 14508 19137 x5+2*x4+3*x3+4*x2+5*x+6ans=8256 10887 14514 191373.3.多项式乘法和除法多项式乘法和除法convconv函函数数进进行行卷卷积积(convolution)(convolution)和和多多项项式乘法运算，式乘法运算，C=conv(A,B)对对向向量量A A和和B B进进行行卷卷积积，返返回回结结果果为为长长度度为为length(A)+length(B)-1length(A)+length(B)-1的向量的向量;如如果果A A和和B B是是多多项项式式的的系系数数向向量量，它它们们的的卷卷积积相相当于两个多项式相乘。当于两个多项式相乘。向量卷积向量卷积所谓两个向量卷积，就是多项式乘法。所谓两个向量卷积，就是多项式乘法。比如:p=1 2 3,q=1 1 是两个向量，p p和和q q的卷积如下：的卷积如下：把p的元素作为一个多项式的系数按降幂排列，写出对应的多项式：x2+2x+3;把q的元素也作为多项式的系数按降幂排列，写出对应的多项式：x+1；对两个多项式相乘，(x2+2x+3)(x+1)=x3+3x+5x2+3取系数，所以p和q卷积的结果就是1 3 5 3。3.3.多项式乘法和除法多项式乘法和除法deconv函数进行反卷积(deconvolution)(deconvolution)和多项式除法运算，l Q,R=deconv(B,A)对向量对向量A A和和B B进行反卷积运算，返回结果为向量进行反卷积运算，返回结果为向量Q Q和余量和余量R;R;如如果果A A和和B B是是多多项项式式的的系系数数向向量量，它它们们的的反反卷卷积积相相当当于于两两个多项式相除，个多项式相除，A A是被除数，是被除数，B B是除数，是除数，Q Q是商，是商，R R是余数。是余数。4.4.多项式的微积分多项式的微积分（1 1）多项式的微分）多项式的微分polyderpolyder函数计算多项式的微分函数计算多项式的微分polyder(P)polyder(P)，返回多项式，返回多项式P P微分的系数向量；微分的系数向量；polyder(A,B)polyder(A,B)，返回多项式，返回多项式A*BA*B微分的系数向量；微分的系数向量；Q,D=polyder(B,A)Q,D=polyder(B,A)，返返回回多多项项式式B/AB/A微微分分的的系系数数向量。向量。实 例如何用如何用MATLABMATLAB对一个已知的函数：对一个已知的函数：y=3*x3+0.5x2+7*x-0.09y=3*x3+0.5x2+7*x-0.09进行求导，并分进行求导，并分别作出求导前和求导后的相应曲线。别作出求导前和求导后的相应曲线。P=3 0.5 7-0.09,X=polyder(P),poly2str(X,x)P=3.0000 0.5000 7.0000 -0.0900X=9 1 7ans=9 x2+x+7（2 2）多项式的积分）多项式的积分polyintpolyint函数计算多项式的不定积分，polyint(P,C)polyint(P,C)P为多项式系数向量；C为不定积分常数项，为标量。polyint(P)polyint(P)，假设 C=0.返回值为多项式不定积分的系数向量。5.5.多项式的根和由根创建多项式多项式的根和由根创建多项式（1）多项式的根roots函数用于求多项式的根，其用法如下：r=roots(C)返回以C向量为系数的多项式的所有根r。实 例求方程求方程 3*x3+0.5x2+7*x-0.09=03*x3+0.5x2+7*x-0.09=0的根。的根。P=3 0.5 7-0.09,x=roots(P)P=3.0000 0.5000 7.0000 -0.0900 x=-0.0898+1.5256i -0.0898-1.5256i 0.0128（2）由根创建多项式poly函数实现由根创建多项式，其用法如下：p=poly(r)，输入变量r是多项式的所有根，返回值为多项式的系数向量，与roots函数是两个可逆的过程；p=poly(A)，输入变量A是方阵，返回值为A的特征多项式的系数向量。部分分式展开部分分式展开又叫部分分式分解，是将有理函数分解成许多次数较低有理函数和的形式，来降低分子或分母多项式的次数，例如：分解后的分式需满足以下条件：分式的分母需为不可约多项式；分式的分子多项式次数需比其分母多项式次数要低。6.6.多项式部分分式展开多项式部分分式展开6.6.多项式部分分式展开多项式部分分式展开residue函数将多项式之比按部分分式展开，其用法如下：lR,P,K R,P,K=residue(B,A)residue(B,A)，求求多多项项式式B/AB/A的的部分分式展开部分分式展开(分解分解)lB,A=residue(R,P,K)B,A=residue(R,P,K)，从从部部分分分分式式得得到到多多项式向量。项式向量。R,P,K=residue(B,A)R,P,K=residue(B,A)B(x)R(1)R(2)R(n)B(x)R(1)R(2)R(n)-=-+-+.+-=-+-+.+-+K(x)+K(x)A(x)x-P(1)x-P(2)x-P(n)A(x)x-P(1)x-P(2)x-P(n)6.6.多项式部分分式展开多项式部分分式展开7.7.多项式曲线拟合多项式曲线拟合l已知在点集上的函数值,构造一个解析函数(其图形为一曲线)使在原离散点上尽可能接近给定的值，这一过程称为曲线拟合。polyfitpolyfit函函数数采采用用最最小小二二乘乘法法对对给给定定数数据据进进行行多多项项式式拟合，其用法如下：拟合，其用法如下：p=polyfit(x,y,n)p=polyfit(x,y,n)p,s=polyfit(x,y,n)p,s=polyfit(x,y,n)说明：说明：x,y为数据点，n为多项式阶数，返回p为幂次从高到低的多项式系数向量p。x必须是单调的。矩阵s用于生成预测值的误差估计。实 例 1给出如下数据的拟合曲线，给出如下数据的拟合曲线，x=0.5,1.0,1.5,2.0,2.5,3.0 x=0.5,1.0,1.5,2.0,2.5,3.0，y=1.75,2.45,3.81,4.80,7.00,8.60y=1.75,2.45,3.81,4.80,7.00,8.60MATLAB程序如下：x=0.5,1.0,1.5,2.0,2.5,3.0;y=1.75,2.45,3.81,4.80,7.00,8.60;p=polyfit(x,y,2)x1=0.5:0.05:3.0;y1=polyval(p,x1);plot(x,y,*r,x1,y1,-b)计算结果为：p=0.5614 0.8287 1.1560即所得多项式为y=0.5614x2+0.08287x+1.155608.8.多项式曲线拟合图形用户接口多项式曲线拟合图形用户接口l曲线拟合的图形用户接口可通过图形窗口的【Tools】菜单中【Basic Fitting】选项启动。x=0:0.1:10;y=0.25*x+20*sin(x);plot(x,y,ro)最小二乘法曲线拟合最小二乘法曲线拟合利用利用lsqcurvefit函数进行任意表达式曲线拟合函数进行任意表达式曲线拟合原理：原理：min sum(fun(x,xdata)-ydata).2用法：用法：x=x=lsqcurvefitlsqcurvefit(fun,x0,xdata,ydata)(fun,x0,xdata,ydata)x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub)x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub)x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)x,resnormx,resnorm=lsqcurvefit(fun,x0,xdata,ydata,.)=lsqcurvefit(fun,x0,xdata,ydata,.)x,resnorm,residual,exitflagx,resnorm,residual,exitflag=lsqcurvefit(fun,x0,xda=lsqcurvefit(fun,x0,xdata,ydata,.)ta,ydata,.)实实 例例xdata=0:0.1:10;ydata=3*sin(xdata)+6;fun=inline(x(1)*sin(xdata)+x(2),x,xdata);x=lsqcurvefit(fun,0 0,xdata,ydata)yfit=feval(fun,x,xdata);plot(xdata,ydata,ro,xdata,yfit)legend(ydata,yfit)5.2 5.2 插插 值值 运运 算算 插值运算是根据数据的分布规律，找到一个函数表达插值运算是根据数据的分布规律，找到一个函数表达式可以连接已知的各点，并用此函数表达式预测两点式可以连接已知的各点，并用此函数表达式预测两点之间任意位置上的函数值。之间任意位置上的函数值。5.2.1 5.2.1 一维插值一维插值5.2.2 5.2.2 二维插值二维插值插值是根据已知输入/输出数据集和当前输入估计输出值。MATLAB提供大量的插值函数，如下表所示。插值函数5.2.1 5.2.1 一维插值一维插值一维插值就是对函数 y=f(x)进行插值，一维插值的原理如下图所示。interp1interp1函数实现一维插值，其用法如下：函数实现一维插值，其用法如下：lyi=interp1(x,y,xi)yi=interp1(x,y,xi)，x x、y y是已知数据集且具有相同长度的是已知数据集且具有相同长度的向量；向量；lyi=interp1(y,xi)yi=interp1(y,xi)，默认，默认x x为为1:n1:n，其中，其中n n为向量为向量y y的长度；的长度；lyi=interp1(x,y,xi,method)yi=interp1(x,y,xi,method)methodmethod用于指定插值的方法，有用于指定插值的方法，有nearest,linear,spline,pchip,cubic,v5cubic等方法。等方法。各种插值方法在速度上各种插值方法在速度上,Nearest,Nearest最快最快,然后是然后是LinearLinear再到再到Cubic,Cubic,最慢的是最慢的是Spline.Spline.但是精度和曲线的平滑度恰好相反但是精度和曲线的平滑度恰好相反,Nearest,Nearest甚至不连续甚至不连续系统默认的是系统默认的是LinearLinear。内插运算与外插运算内插运算与外插运算(1)只对已知数据点集内部的点进行的插值运算称为内插，可比较准确的估测插值点上的函数值；(2)当插值点落在已知数据集的外部时的插值称为外插，要估计外插函数值很难。y=interp1(x,y,xi,method,extrap,extrapval)代表可选。MATLAB对已知数据集外部点上函数值的预测都返回NaN，但可通过为interp1函数添加extrap参数指明用于外插，MATLAB的外插结果偏差较大。实 例Sw Kro Krw0.290 1 0.000 0.356 0.63 0.017 0.449 0.35 0.0690.496 0.27 0.1130.543 0.2 0.1650.590 0.133 0.2270.684 0.05 0.3720.730 0.01 0.4560.780 0 0.546问题：问题：已知相渗曲线（如图所示），原油的密度为0.9，粘度为9，地层水的密度为1，粘度为0.45，计算一定饱和度时的含水率。实 例load xiangshen.datSwg=xiangshen(:,1);Krog=xiangshen(:,2);Krwg=xiangshen(:,3);plot(Swg,Krog,r-o,Swg,Krwg,b-o)legend(Kro,Krw)miuo=9;rouo=0.9;miuw=0.45;rouw=1;load waterSw=WAT;Kro=interp1(Swg,Krog,Sw);Krw=interp1(Swg,Krwg,Sw);fw=Krw*rouw/miuw./(Krw*rouw/miuw+Kro*rouo/miuo);figure,surf(Sw),colorbar,view(2),axis equal,axis tightfigure,surf(fw),colorbar,view(2),axis equal,axis tight5.2.2 5.2.2 二维插值二维插值二二维维插插值值是是对对两两变变量量的的函函数数z=f(x,y)z=f(x,y)进进行行插插值值，二二维插值的原理如下图所示：维插值的原理如下图所示：interp2函数实现二维插值,其用法如下：lzi=interp2(x,y,z,xi,yi),x,y,zzi=interp2(x,y,z,xi,yi),x,y,z为为原原始始数数据据，返返回回值值zizi是是插插值结果；值结果；lzi=interp2(z,xi,yi),zi=interp2(z,xi,yi),若若z=nz=nm,m,则则x=1:n,y=1:mx=1:n,y=1:m；lzi=interp2(x,y,z,xi,yi,method),zi=interp2(x,y,z,xi,yi,method),method用于指定插值的方法,外插使用extrapval 。已知数据点集(x x,y y)必须是栅格格式，一般用meshgrid函数产生。interp2要求(x x,y y)必须是严格单调的并且是等间距的，如果(x x,y y)不是等间距的，会将且变换为等间距形式，如果已知是等间距的，可在method参数前加星号，如果：*cubic。x,y,z=peaks(10);xi,yi=meshgrid(-3:.1:3,-3:.1:3);zi=interp2(x,y,z,xi,yi);subplot(1,2,1);mesh(xi,yi,zi);subplot(1,2,2);mesh(x,y,z);实实 例例5.3 数 据 分 析 5.3.1 5.3.1 基本数据分析函数基本数据分析函数 5.3.2 5.3.2 协方差和相关系数矩阵协方差和相关系数矩阵 5.3.3 5.3.3 有限差分和梯度有限差分和梯度 5.3.4 信号滤波和卷积 5.3.5 傅立叶变换 MATLAB提供大量数据分析的函数，首先给出如下约定：l一维数据分析时，数据可以用行向量或者列向量来表示；l二维数据分析时，数据可以用多个向量或者二维矩阵来表示。5.3.1 基本数据分析函数 1 1最大值、最小值、平均值、最大值、最小值、平均值、中间值、元素求和中间值、元素求和 2 2标准差和方差标准差和方差 3 3元素排序元素排序基本数据分析函数基本数据分析函数 续表1.最值、平均值、中间值、和 运行结果如下，并如下图所示。运行结果如下，并如下图所示。2.标准差和方差3.3.元素排序元素排序MATLAB提供对实数、复数和字符串的排序函数。lsortsort函数实现对数值的升序排列；函数实现对数值的升序排列；lsortrowssortrows函数实现对行的升序排序。函数实现对行的升序排序。实实 例例 X=magic(4)X=magic(4)X=16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1 Y,I=sort(X)Y,I=sort(X)Y=4 2 3 1 5 7 6 8 9 11 10 12 16 14 15 13I=4 1 1 4 2 3 3 2 3 2 2 3 1 4 4 1 结果解释：对矩阵X的每列按升序排列得到新矩阵Y，并返回序数矩阵。X=magic(4)X=magic(4)X=16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1 Y,I=sort(X,2)Y,I=sort(X,2)Y=2 3 13 16 5 8 10 11 6 7 9 12 1 4 14 15I=2 3 4 1 1 4 3 2 3 2 1 4 4 1 2 3结果解释：对矩阵X的每行按升序排列得到新矩阵Y，并返回序数矩阵。5.3.2 5.3.2 协方差和相关系数矩阵协方差和相关系数矩阵 cov函数计算随机变量的协方差矩阵，其用法如下：lC=cov(X)，计算X代表的随机变量的协方差矩阵；lC=cov(X,Y)，X和Y必须是具有相同长度的向量，等同于cov(X(:)Y(:).corrcoef函数计算随机变量的相关系数矩阵，其用法如下：lR=corrcoef(X)，返回X代表的随机变量的相关系数矩阵；lR=corrcoef(X,Y)，X和Y必须是具有相同长度的列向量，等同于R=corrcoef(X Y).5.3.3 5.3.3 有限差分和梯度有限差分和梯度diff函数计算差分，其用法如下：lY=diff(X)，X可以是向量或矩阵，如果X是向量，Y=X(2)-X(1)X(3)-X(2).X(n)-X(n-1)；如果X是矩阵，Y=X(2:n,:)-X(1:n-1,:)。lY=diff(X,n)，计算n阶差分；lY=diff(X,n,dim)，计算在dim维上的n阶差分。实实 例例问题：已知压力计算压力导数 load pwdvstd loglog(tD,pwD)hold on loglog(tD(1:end-1),tD(1:end-1).*diff(pwD)./diff(tD)gradient函数计算梯度，其用法如下：lFX=gradient(F)，返回F在x方向上的梯度；lFX,FY=gradient(F)，FX是F在x方向的近似偏导数，FY是F在y方向的近似偏导数；lFx,Fy,Fz,.=gradient(F)，返回N个方向的近似偏导数；l.=gradient(F,h)，h用于指定所有方向上自变量的间距；l.=gradient(F,h1,h2,.)，用多个标量来指定各个方向上自变量的间距。load prsDX=50;DY=50;NX=101;NY=101;X,Y=meshgrid(DX/2:DX:DX*101-DX/2,DY/2:DY:DY*101-DY/2);K=0.8;miu=0.4;U,V=gradient(pressure,DX,DY);U=-K/miu*U;V=-K/miu*V;starty1=50*DY:2.5:51*DY;startx1=50*DX*ones(size(starty1);starty2=50*DY:2.5:51*DY;startx2=51*DX*ones(size(starty2);startx3=50*DX:2.5:51*DX;starty3=50*DY*ones(size(startx3);startx4=50*DX:2.5:51*DX;starty4=51*DY*ones(size(startx4);startx=startx1 startx2 startx3 startx4;starty=starty1 starty2 starty3 starty4;streamline(X,Y,U,V,startx,starty)axis equalaxis tightbox on实实 例例5.3.4 信号滤波和卷积1一维数字滤波 2信号卷积 3去除信号直流或线性成分 MATLAB提供如下表所示的信号滤波和卷积的函数。信号滤波和卷积函数 1一维数字滤波 函数filter()实现一维数字滤波，该函数的具体使用方法如下：y=filter(b,a,X)，X为用于滤波的数据，Y为数据X通过滤波器之后的值；y,zf=filter(b,a,X)，附加返回一个表示数据延迟时间的量zf；y,zf=filter(b,a,X,zi)，zi为初始数据延迟，zf等于最终数据延迟；y=filter(b,a,X,zi,dim)，在dim维上进行数据滤波。2信号卷积 函数conv计算卷积。运行结果如下图所示。3去除信号直流或线性成分 detrend函数实现去除信号中的直流或者线性成分，其具体使用方法如下：y=detrend(x)，如果x是一个向量，从信号x中减去线性成分；如果x是一个矩阵，去除x所有列中的线性成分；y=detrend(x,constant)，如果x是一个向量，减去信号中的直流成分；如果x是一个矩阵，去除所有列中的直流成分；y=detrend(x,linear,bp)，从信号x中减去分段线性函数。5.3.5 傅立叶变换1一维傅立叶变换和逆变换2二维傅立叶变换和逆变换 傅立叶变换既可以对连续信号进行变换，也可以对离散信号进行变换。本小节只介绍离散傅立叶变换。傅立叶变换函数 1一维傅立叶变换和逆变换fft函数实现一维离散傅立叶变换，其用法如下：Y=fft(X)，如果X是向量，返回向量X的傅立叶变换；如果X是矩阵，函数对矩阵X的每一列进行傅立叶变换；Y=fft(X,n)，用输入n指定傅立叶变换的长度；Y=fft(X,dim)，在dim维上进行傅立叶变换；Y=fft(X,n,dim)，在dim维上进行傅立叶变换，并指定傅立叶变换的长度。ifft函数实现一维离散傅立叶逆变换，其用法与fft函数类似，只是添加一个选项。y=ifft(.,symmetric)；y=ifft(.,nonsymmetric)。2二维傅立叶变换和逆变换 fft2函数实现二维傅立叶变换，用ifft2函数来实现二维傅立叶逆变换。fft2函数用法如下：Y=fft2(X)，X是矩阵，对矩阵X进行二维傅立叶变换；Y=fft2(X,m,n)，m和n指定傅立叶变换的长度。5.4 功 能 函 数 1函数的表示 2函数画图 3函数最小值和零点 4数值积分 5在功能函数中使用含参函数 函数可以通过以下方式来表示：M文件；匿名函数；函数inline()。1.函数的表示2函数画图 MATLAB提供函数画图的函数如下表所示。函数画图的函数 fplot函数为例介绍画图函数的用法，其用法如下：fplot(function,limits)，function为待画图的函数，limits是横坐标数值范围或横纵坐标数值范围；fplot(function,limits,LineSpec)，LineSpec指定画图的线条属性；fplot(function,limits,tol)，tol指定画图相对精度；fplot(function,limits,tol,LineSpec)，指定画图的线条属性和画图相对精度。运行结果如下图所示。3函数最小值和零点求函数的最小值和零点的函数，如下表所示。求函数最小值和零点（1）求一元函数最小值fminbnd函数求一元函数在给定区间内的最小值，其具体使用方法如下：x=fminbnd(fun,x1,x2)，在区间x1 x2内寻找函数最小值；x=fminbnd(fun,x1,x2,options)，使 用options选项来指定的优化器的参数；x,fval=fminbnd(.)，附加返回函数最小值。（2）求多元函数的最小值fminsearch函数求多元函数的最小值。其用法如下：x=fminsearch(fun,x0)，在初始x0附近寻找局部最小值；x=fminsearch(fun,x0,options)，使 用options选项来指定优化器的参数；x,fval=fminsearch(.)，附加返回函数最小值。（3）求一元函数的零点fzero函数求一元函数的零点，其用法如下：x=fzero(fun,x0)，在x0点附近寻找函数的零点；x=fzero(fun,x0,x1)，在x0,x1区间内寻找函数的零点；x=fzero(fun,x0,options)，用options指定寻找零点的优化器参数；x,fval=fzero(.)，附加自变量为x时的函数值。（4）优化器参数optimset函数设定优化器参数，其用法如下：options=optimset(param1,value1,param2,value2,.)，用参数名和对应的参数值设定优化器的参数；optimset，显示优化器的所有参数名和有效的参数值；options=optimset，返回一个优化器的结构体，所有的域的值都为；options=optimset(optimfun)，返回函数optimfun对应的优化器参数。options=optimset(oldopts,param1,value1,.)，在原优化器参数oldopts的基础上，改动指定优化器参数；options=optimset(oldopts,newopts)，用newopts的所有非空参数覆盖oldopts中的值。在optimset函数中常用的优化器参数如下表所示。优化器参数 optimget函数得到目前优化器的参数，其具体使用方法如下：val=optimget(options,param)，返回优化器参数param的值；val=optimget(options,param,default)，返回优化器参数param的值。4数值积分 MATLAB提供一些的数值积分函数，如下表所示。数值积分函数（1）一元函数的数值积分quad函数和quadl函数用来计算一元函数的积分。quad函数的用法如下：q=quad(fun,a,b)，计算函数fun在a b区间内的定积分；q=quad(fun,a,b,tol)，以绝对误差容限tol计算函数fun在a b区间内的定积分；q=quad(fun,a,b,tol,trace)，当trace为非零值时，显示迭代过程的中间值。（2）矢量数值积分 矢量数值积分等价于多个一元定积分。（3）二重和三重积分dblquad函数计算二重积分。其具体使用方法如下：q=dblquad(fun,xmin,xmax,ymin,ymax)，计算二元函数的二重积分；q=dblquad(fun,xmin,xmax,ymin,ymax,tol)，用tol指定绝对计算精度；q=dblquad(fun,xmin,xmax,ymin,ymax,tol,method)，用method指定计算一维积分时采用的函数。5在功能函数中使用含参函数 功能函数中的含参函数的两种解决方法：嵌套函数匿名函数 编写M文件的函数时：首先将含参函数的参数作为输入；其次在其中调用功能函数，形成嵌套；最后通过调用该函数进行计算。（1）用嵌套函数提供函数参数运行结果如下图所示。（2）用匿名函数提供函数参数 用匿名函数提供函数参数的具体步骤如下：创建一个含参函数，并保存为M文件格式；调用功能函数的M文件中给参数赋值；用含参函数创建匿名函数；把匿名函数句柄传递给功能函数计算。5.5 微分方程组数值解 5.5.1 常微分方程组的初值问题 5.5.2 延迟微分方程的问题 5.5.3 常微分方程组的边界问题 在MATLAB中，可以计算微分方程数值解，如：常微分方程组的初值问题；延迟微分方程的问题；常微分方程组的边界问题。5.5.1 常微分方程组的初值问题 1显式常微分方程组 2设置解法器参数 3线性隐式常微分方程组 4完全隐式常微分方程组在MATLAB中可以计算以下初值问题的数值解。显式常微分方程组；线性隐式常微分方程组；完全隐式常微分方程组。1显式常微分方程组在MATLAB中，用函数实现不同的解法，如下表所示。常微分方程组解法对比 2设置解法器参数odeset函数设定解法器参数，其用法如下：options=odeset(name1,value1,name2,value2,)，用参数名和相应参数值设定解法器的参数；options=odeset(oldopts,name1,value1,)，修改原来的解法器options结构体oldopts；options=odeset(oldopts,newopts)，合并两个解法器options结构体oldopts和newopts；odeset，显示所有的参数值和它们的默认值。常微分方程组解法器参数3线性隐式常微分方程组 线性隐式常微分方程组可以利用解法器参数options来求解。运行结果如下图所示。4完全隐式常微分方程组ode15i函数求解完全隐式常微分方程组，其用法如下：t,Y=ode15i(odefun,tspan,y0,yp0)；t,Y=ode15i(odefun,tspan,y0,yp0,options)。decic函数得到自洽初始值，其用法如下：y0mod,yp0mod=decic(odefun,t0,y0,fixed_y0,yp0,fixed_yp0)；y0mod,yp0mod=decic(odefun,t0,y0,fixed_y0,yp0,fixed_yp0,options)。运行结果如下图所示。6.5.2 延迟微分方程的问题dde23函数求解延迟微分方程组，其用法如下：sol=dde23(ddefun,lags,history,tspan)；sol=ode23(ddefun,lags,history,tspan,option)，option结构体用于设置解法器的参数。dde23函数的返回值是一个结构体，它包含7个属性，其中重要的5个属性如下：sol.x，dde23选择计算的时间点；sol.y，在时间点x上的解y(x)；sol.yp，在时间点x上解的一阶导数y(x)；sol.history，方程初始值；sol.solver，解法器的名字dde23；若需得到tint时刻的解，可以使用函数deval，即yint=deval(sol,tint)。5.5.3 常微分方程组的边界问题bvp4c 函数的用法如下：sol=bvp4c(odefun,bcfun,solinit)，odefun代表常微分方程组的函数，bcfun是描述边界条件的函数，solinit是对方程解的猜测解；sol=bvp4c(odefun,bcfun,solinit,options)，使用options结构体来设定解法器的参数。运行结果如下图所示。习 题