第5章数值计算.

第五章数值计算目录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([102],'s')ans=s^2+22.多项式求值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)P=123456ans=x^5+2x^4+3x^3+4x^2+5x+6x=5ans=4881例子：P=1:6,poly2str(P,'x'),x=[23;45],polyval(P,x)P=123456ans=x^5+2x^4+3x^3+4x^2+5x+6x=2345ans=12054318184881Y=polyval(P,X),把矩阵或向量X中的每个元素逐个代入多项式中进行计算；Y=polyvalm(P,X)，把矩阵X作为整体代入多项式中进行计算，X必须为方阵。例子：P=1:6,poly2str(P,'x'),x=[23;45],polyval(P,x)P=123456ans=x^5+2x^4+3x^3+4x^2+5x+6x=2345ans=12054318184881例子：P=1:6,poly2str(P,'x'),x=[23;45],polyvalm(P,X)P=123456ans=x^5+2x^4+3x^3+4x^2+5x+6x=2345ans=8256108811450819137x^5+2*x^4+3*x^3+4*x^2+5*x+6ans=82561088714514191373.多项式乘法和除法conv函数进行卷积(convolution)和多项式乘法运算，C=conv(A,B)对向量A和B进行卷积，返回结果为长度为length(A)+length(B)-1的向量;如果A和B是多项式的系数向量，它们的卷积相当于两个多项式相乘。向量卷积所谓两个向量卷积，就是多项式乘法。比如:p=[123],q=[11]是两个向量，p和q的卷积如下：•把p的元素作为一个多项式的系数按降幂排列，写出对应的多项式：x^2+2x+3;•把q的元素也作为多项式的系数按降幂排列，写出对应的多项式：x+1；•对两个多项式相乘，(x^2+2x+3)×(x+1)=x^3+3x+5x^2+3•取系数，所以p和q卷积的结果就是[1353]。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微分的系数向量。实例•如何用MATLAB对一个已知的函数：y=3*x^3+0.5x^2+7*x-0.09进行求导，并分别作出求导前和求导后的相应曲线。P=[30.57-0.09],X=polyder(P),poly2str(X,'x')P=3.00000.50007.0000-0.0900X=917ans=9x^2+x+7（2）多项式的积分polyint函数计算多项式的不定积分，polyint(P,C)P为多项式系数向量；C为不定积分常数项，为标量。polyint(P)，假设C=0.返回值为多项式不定积分的系数向量。5.多项式的根和由根创建多项式（1）多项式的根roots函数用于求多项式的根，其用法如下：r=roots(C)返回以C向量为系数的多项式的所有根r。实例•求方程3*x^3+0.5x^2+7*x-0.09=0的根。P=[30.57-0.09],x=roots(P)P=3.00000.50007.0000-0.0900x=-0.0898+1.5256i-0.0898-1.5256i0.0128（2）由根创建多项式poly函数实现由根创建多项式，其用法如下：p=poly(r)，输入变量r是多项式的所有根，返回值为多项式的系数向量，与roots函数是两个可逆的过程；p=poly(A)，输入变量A是方阵，返回值为A的特征多项式的系数向量。•部分分式展开又叫部分分式分解，是将有理函数分解成许多次数较低有理函数和的形式，来降低分子或分母多项式的次数，•例如：•分解后的分式需满足以下条件：–分式的分母需为不可约多项式；–分式的分子多项式次数需比其分母多项式次数要低。6.多项式部分分式展开6.多项式部分分式展开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)6.多项式部分分式展开7.多项式曲线拟合已知在点集上的函数值,构造一个解析函数(其图形为一曲线)使在原离散点上尽可能接近给定的值，这一过程称为曲线拟合。polyfit函数采用最小二乘法对给定数据进行多项式拟合，其用法如下：•p=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]，y=[1.75,2.45,3.81,4.80,7.00,8.60]MATLAB程序如下：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.56140.82871.1560即所得多项式为y=0.5614x^2+0.08287x+1.155608.多项式曲线拟合图形用户接口曲线拟合的图形用户接口可通过图形窗口的【Tools】菜单中【BasicFitting】选项启动。x=0:0.1:10;y=0.25*x+20*sin(x);plot(x,y,'ro')最小二乘法曲线拟合•利用lsqcurvefit函数进行任意表达式曲线拟合•原理：minsum{(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;ydata=3*sin(xdata)+6;fun=inline('x(1)*sin(xdata)+x(2)','x','xdata');x=lsqcurvefit(fun,[00],xdata,ydata)yfit=feval(fun,x,xdata);plot(xdata,ydata,'ro',xdata,yfit)legend('ydata','yfit')5.2插值运算插值运算是根据数据的分布规律，找到一个函数表达式可以连接已知的各点，并用此函数表达式预测两点之间任意位置上的函数值。5.2.1一维插值5.2.2二维插值插值是根据已知输入/输出数据集和当前输入估计输出值。MATLAB提供大量的插值函数，如下表所示。插值函数5.2.1一维插值一维插值就是对函数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’等方法。各种插值方法在速度上,Nearest最快,然后是Linear再到Cubic,最慢的是Spline.但是精度和曲线的平滑度恰好相反,Nearest甚至不连续~~系统默认的是Linear。内插运算与外插运算(1)只对已知数据点集内部的点进行的插值运算称为内插，可比较准确的估测插值点上的函数值；(2)当插值点落在已知数据集的外部时的插值称为外插，要估计外插函数值很难。y=interp1([x,]y,xi,[method],['extrap'],[extrapval])[]代表可选。•MATLAB对已知数据集外部点上函数值的预测都返回NaN，但可通过为interp1函数添加'extrap'参数指明用于外插，MATLAB的外插结果偏差较大。实例SwKroKrw0.29010.0000.3560.630.0170.4490.350.0690.4960.270.1130.5430.20.1650.5900.1330.2270.6840.050.3720.7300.010.4560.78000.546问题：已知相渗曲线（如图所示），原油的密度为0.9，粘度为9，地层水的密度为1，粘度为0.45，计算一定饱和度时的含水率。实例loadxiangshen.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;loadwaterSw=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),axisequal,axistightfigure,surf(fw),colorbar,view(2),axisequal,axistight5.2.2二维插值二维插值是对两变量的函数z=f(x,y)进行插值，二维插值的原理如下图所示：interp2函数实现二维插值,其用法如下：zi=interp2(x,y,z,xi,yi),x,y,z为原始数据，返回值zi是插值结果；zi=interp2(z,xi,yi),若z=n×m,则x=1: