Matlab编程技术：第5章 数值计算

1、第五章 数值计算,吴明录 2011-3-26,目 录,5.1 多项式运算 5.2 插 值 运 算 5.3 数 据 分 析 5.4 功 能 函 数 5.5 微分方程组数值解 习 题,5.1 多项式运算,1.多项式Matlab表示法 2.多项式求值 3.多项式乘法和除法 4.多项式的微积分 5.多项式的根和由根创建多项式 6.多项式部分分式展开 7.多项式曲线拟合 8.多项式曲线拟合图形用户接口,Matlab提供了下列关于多项式的函数,1.多项式表示法,MATLAB采用行向量表示多项式系数，多项式系数按降幂排列。 poly2str函数将多项式系数向量转换为完整形式。 poly2str(1 0 2,

2、s) ans = s2 + 2,2.多项式求值,polyval函数计算多项式的值，其用法为： Y = polyval(P,X), P多项式系数行向量 X代入多项式的值 其意义为： 当P为长度为n+1的行向量 p1,p2, .,pn,pn+1 时， Y = p1xn + p2xn-1 + . + pnx + pn+1,例子： P=1:6 poly2str(P,x) x=5, polyval(P,x) A=2 3;4 5 polyval(P,A,Y = polyval(P,X)： 把矩阵或向量X中的每个元素逐个代入多项式中进行计算； Y = polyvalm(P,X)： 把矩阵X作为整体代入多项式

3、中进行计算，X必须为方阵,例子： P=1:6, poly2str(P,X), X=2 3;4 5,polyval(P,X) polyvalm(P,X,比较,3.多项式乘法和除法,conv函数进行卷积(convolution)和多项式乘法运算 C = conv(A, B) A，B多项式的系数向量 C长度为length(A)+length(B)-1的向量,向量卷积,向量卷积就是多项式乘法 对于 p=1 2 3, q=1 1 把p的元素作为一个多项式的系数按降幂排列，写出对应的多项式：x2+2x+3; 把q的元素也作为多项式的系数按降幂排列，写出对应的多项式：x+1； 对两个多项式相乘， (x2+2

4、x+3)(x+1)=x3+3x+5x2+3 取系数，p和q卷积的结果就是1 3 5 3,3.多项式乘法和除法,deconv函数进行反卷积(deconvolution)和多项式除法运算 Q,R = deconv(B,A) A，B多项式向量； Q结果向量； R余量。 A和B的反卷积相当于两个多项式相除，A是被除数，B是除数，Q是商，R是余数,4.多项式的微积分,1）多项式的微分polyder polyder(P) 返回多项式P微分的系数向量； polyder(A,B) 返回多项式A*B微分的系数向量； Q,D=polyder(B,A) 返回多项式B/A微分的系数向量Q/D,2）多项式的不定积分po

5、lyint,polyint(P,C) P多项式系数向量； C不定积分常数项，标量。 polyint(P) C=0 返回值为多项式不定积分的系数向量,实 例,如何对一个已知的函数：y=3*x3+0.5x2+7*x-0.09进行求导，并分别作出求导前后在-10 10的曲线。 P=3 0.5 7 -0.09 poly2str(P,x) dP=polyder(P) poly2str(dP,x) x=-10:0.5:10； plot(x,polyval(P,x),r-,x, polyval(dP,x),b,5.多项式的根和由根创建多项式,1）多项式的根roots r = roots(C) 返回以C向量为

6、系数的多项式的所有根r,例子：求方程 3*x3+0.5x2+7*x-0.09=0的根。 P=3 0.5 7 -0.09 x=roots(P,2）由根创建多项式poly 与roots函数是两个可逆的过程 p=poly(r) r多项式的所有根 p多项式的系数向量 p=poly(A) A方阵 返回值为A的特征多项式的系数向量,又叫部分分式分解，是将有理函数分解成许多次数较低有理函数和的形式，来降低分子或分母多项式的次数 例如： 分解后的分式需满足以下条件： 分式的分母需为不可约多项式； 分式的分子多项式次数需比其分母多项式次数要低,6.多项式部分分式展开,6.多项式部分分式展开,将多项式之比按部分分

7、式展开residue R,P,K = residue(B,A) 求多项式B/A的部分分式展开(分解) B,A=residue(R,P,K) 从部分分式得到多项式向量,R,P,K = residue(B,A) B(x) R(1) R(2) R(n) - = - + - + . + - + K(x) A(x) x - P(1) x - P(2) x - P(n,7.多项式曲线拟合,曲线拟合已知在点集上的函数值,构造一个解析函数(其图形为一曲线)使在原离散点上尽可能接近给定的值的过程。 polyfit函数采用最小二乘法对给定数据进行多项式拟合 p=polyfit(x,y,n) p,s= polyfi

8、t(x,y,n) x,y数据点， x必须是单调的 n多项式阶数 p多项式系数向量 s矩阵，用于生成预测值的误差估计,实 例,给出如下数据的拟合曲线 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 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,8.多项式曲线拟合图形用户界面,最小二乘法曲线拟合,利用lsqcur

9、vefit函数进行任意表达式曲线拟合 原理： min sum (fun(x,xdata)-ydata).2 用法： x=lsqcurvefit(fun,x0,xdata,ydata) x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub) x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options) x,resnorm=lsqcurvefit(fun,x0,xdata,ydata,.) x,resnorm,residual,exitflag=lsqcurvefit(fun,x0,xdata,ydata,.,实 例,xdata=0:0.1

10、:10; ydata=3*sin(xdata)+6; fun=inline(x(1)*sin(xdata)+x(2),x,xdata); % fun=(x,xdata) x(1)*sin(xdata)+x(2); 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.1 一维插值 5.2.2 二维插值,插值函数,5.2.1 一维插值,一维插

11、值就是对函数 y=f(x) 进行插值,一维插值interp1 yi=interp1(x,y,xi) x、y已知数据集且具有相同长度的向量 yi=interp1(y,xi) 默认x为1:n，其中n为向量y的长度； yi=interp1(x,y,xi,method) method用于指定插值的方法，有nearest, linear,spline,pchip, cubic,v5cubic等方法,内插运算与外插运算,1)只对已知数据点集内部的点进行的插值运算称为内插，可比较准确的估测插值点上的函数值； (2)当插值点落在已知数据集的外部时的插值称为外插，要估计外插函数值很难。 y=interp1(x,

12、y,xi,method,extrap,extrapval) 可选。 MATLAB对已知数据集外部点上函数值的预测都返回NaN，但可通过为interp1函数添加extrap参数指明用于外插，MATLAB的外插结果偏差较大,实 例,Sw Kro Krw 0.290 1 0.000 0.356 0.63 0.017 0.449 0.35 0.069 0.496 0.27 0.113 0.543 0.2 0.165 0.590 0.133 0.227 0.684 0.05 0.372 0.730 0.01 0.456 0.780 0 0.546,问题：已知相渗曲线（如图），原油的密度为0.9，粘度为9

13、，地层水的密度为1，粘度为0.45，计算一定饱和度时的含水率,5.2.2 二维插值,二维插值是对两变量的函数z=f(x,y)进行插值,二维插值interp2 zi=interp2(x,y,z,xi,yi) x,y,z原始数据，必须是栅格格式，一般用meshgrid函数产生，(x,y)必须是严格单调的并且是等间距的，如果(x, y)不是等间距的，会将且变换为等间距形式，如果已知是等间距的，可在method参数前加星号，如果：*cubic。 zi插值结果 zi=interp2(z,xi,yi) zi=interp2(x,y,z,xi,yi,method,x,y,z = peaks(10); xi,

14、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.2 协方差和相关系数矩阵 5.3.3 有限差分和梯度 5.3.4 信号滤波和卷积 5.3.5 傅立叶变换,5.3.1 基本数据分析函数,1最大(小)值、平均值、中间值、元素求和 2标准差和方差 3元素排序,基本数据分析函数,续表,1.最值、平均值、中间值、和,运行结果如下,2.标准差和方差,2.标

15、准差（std）和方差（var,向量的标准差有两种定义,s = std(X) 按公式1计算向量X的标准差 s = std(X,flag) 若flag=0按公式1计算标准差；若flag=1则按公式2计算 s = std(X,flag,dim) 沿维数计算标准差，dim=1则沿行计算，dim=2则沿列计算,3.元素排序,MATLAB提供对实数、复数和字符串的排序函数。 sort函数实现对数值的升序排列； sortrows函数实现对行的升序排序,实 例,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 = 4 2 3

16、 1 5 7 6 8 9 11 10 12 16 14 15 13 I = 4 1 1 4 2 3 3 2 3 2 2 3 1 4 4 1 每列按升序排列得到新矩阵Y，并返回序数矩阵,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 = 2 3 13 16 5 8 10 11 6 7 9 12 1 4 14 15 I = 2 3 4 1 1 4 3 2 3 2 1 4 4 1 2 3 每行按升序排列得到新矩阵Y，并返回序数矩阵,5.3.2 协方差和相关系数矩阵,cov函数计算随机变量的协方差矩阵 C=co

17、v(X) 计算X代表的随机变量的协方差矩阵 C=cov(X,Y)等同于cov(X(:) Y(:) X、Y具有相同长度的向量,corrcoef函数计算随机变量的相关系数矩阵 R=corrcoef(X) 返回X代表的随机变量的相关系数矩阵； R=corrcoef(X,Y)等同于R=corrcoef (X Y). X、Y具有相同长度的列向量,5.3.3 有限差分和梯度,diff函数计算差分 Y=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,:)。 Y=diff(X,n) ，

18、计算n阶差分； Y=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函数计算梯度 FX=gradient(F) 返回F在x方向上的梯度 FX,FY=gradient(F) FX、 FYF在x、y方向的近似偏导数 Fx,Fy,Fz,.=gradient(F) 返回N个方向的近似偏导数 .=gradient(F,h) h用于指定所有方向上自变量的间距 .=gradi

19、ent(F,h1,h2,.) 用多个标量来指定各个方向上自变量的间距,load prs DX=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*on

20、es(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 equal axis tight box on,实 例,5.4 功 能 函 数,1函数的表

21、示 2画图函数 3函数最小值和零点 4数值积分,M文件（function） 匿名函数（） 内联函数（inline,1.函数的表示,2. 画图函数,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,opti

22、ons) 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时的函数值,optimget函数得到目前优化器的参数 val = optimget(options,param) 返回优化器参数param的值 val=optimget(options,pa

23、ram,default) 返回优化器参数param的值,4）优化器参数设置,options=optimset(param1,value1, param2, value2,.)用参数名和对应的参数值设定优化器的参数 optimset，显示优化器的所有参数名和有效的参数值； options=optimset，返回一个优化器的结构体，所有的域的值都为； options = optimset(optimfun)，返回函数optimfun对应的优化器参数 options=optimset(oldopts,param1,value1,.)，在原优化器参数oldopts的基础上，改动指定优化器参数； opt

24、ions = optimset(oldopts,newopts)，用newopts的所有非空参数覆盖oldopts中的值,optimset函数可以设置优化器的参数,在optimset函数中常用的优化器参数,优化器参数,4数值积分,1）一元函数的数值积分quad、 quadl,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）矢量数值积分quadv,矢量数值积分等价于多个一元定积

25、分,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.5 微分方程（组）数值解,5.5.1 常微分方程（组）的初值问题 5.5.2 延迟微分方程（组）的问题 5.5.3 常微分方程（组）的边值问题,5.5.1 常微分方程组的初值问题,1显式常微分方程组 2线性隐式常微分方程组 3完

26、全隐式常微分方程组,y=f(x,y,F(x,y,y)=0,1.显式常微分方程组,van der Pol equations,刚性问题,非刚性问题,非刚性问题,function dy = rigid(t,y) dy = zeros(3,1); dy(1) = y(2) * y(3); dy(2) = -y(1) * y(3); dy(3) = -0.51 * y(1) * y(2,options = odeset(RelTol,1e-4,AbsTol,1e-4 1e-4 1e-5); T,Y = ode45(rigid,0 12,0 1 1,options,plot(T,Y(:,1),-,T,Y(:,2),-.,T,Y(:,3),.,范德波耳(van der Pol,刚性问题,function dy = vdp1000(t,y) dy = zeros(2,1); dy(1) =