1、Matlab数值计算第10讲6.1 多项式(polynomial)Application of Matlab Language2013/9/102多项式的MATLAB表达多项式由一个行向量表示该向量元素是该多项式的系数且按降幂次序排列如：多项式x412x325x116由行向量：p=1-12 0 25 116表示。注意，必须包括具有零系数的项。求解多项式的根？roots指令6.1 多项式(polynomial)Application of Matlab Language2013/9/103举例：求解多项式x412x325x116的根p=1-12 025116p=1-12025116r=roots
2、(p)r=11.74732.7028-1.2251+1.4672i-1.2251-1.4672iMATLAB按惯例规定，多项式是行向量，根是列向量6.1 多项式(polynomial)Application of Matlab Language2013/9/104已知多项式的根，求解多项式？能！使用poly指令举例：由上例所得的根求其多项式 pp=poly(r)pp=1.0000-12.0000-0.000025.0000即：x412x325x116116.00006.1 多项式(polynomial)Application of Matlab Language2013/9/105 多项式的乘
3、法(conv指令)举例：多项式a(x)=x32x23x4 和b(x)=x34x29x16的乘积。a=1 2 3 4;b=1 c=conv(a,b)c=162050754 9 16;8464 两个以上的多项式的乘法需要重复使用conv.6.1 多项式(polynomial)Application of Matlab Language2013/9/106多项式加法:MATLAB没有提供进行加法运算的函数。如果两个多项式向量大小相同，标准的数组加法有效。把多项式a(x)与上面给出的b(x)相加。d=a+bd=261220结果:d(x)=2x36x212x20当两个多项式阶次不同，低阶的多项式必须用首
4、零填补，使其 与高阶多项式有同样的阶次。考虑上面多项式c和d相加：e=c+0 0 0 de=162052819684结果:d(x)=x6+6x5+20 x4+52x3+81x2+96x+846.1 多项式(polynomial)Application of Matlab Language2013/9/107 问题：编写一个多项式加法运算的函数文件 明确需求 定义一个函数，需要有两个输入参数（比如：p1、p2），一个输出参数p_out。p1、p2表示两个待计算的多项式。p_out表示两个多项式的求和结果 函数的内部处理 如p1、p2两参数大小相等，则直接相加:p_out=p1+p2 如p1、p2
5、两参数大小不等：If length(p1)length(p2)P2前面要补0元素，使p1、p2两参数大小相等否则P1前面要补0元素，使p1、p2两参数大小相等6.1 多项式(polynomial)Application of Matlab Language2013/9/108函数文件的编写function p_out=poly_sum(p1,p2)%calculate the sum of two polynomialsif length(p1)=length(p2)p_out=p1+p2;elseif length(p1)length(p2)p2=zeros(1,length(p1)-len
6、gth(p2),p2;p_out=p1+p2;elsep1=zeros(1,length(p2)-length(p1),p1;p_out=p1+p2;end6.1 多项式(polynomial)Application of Matlab Language2013/9/109多项式的除法(deconv)举例说明：c(x)=x66x520 x450 x375x284x64 除以b(x)=x34x29x16 c=1 b=1 462050758464;9 16;q,r=deconv(c,b)q=1234r=00000006.1 多项式(polynomial)Application of Matlab
7、Language2013/9/1010多项式的导数(polyder)举例：求b(x)=x34x29x16的导数。b=149 d=polyder(b)d=38916;结果为：3x28x9另外两种形式为：p=polyder(P,Q)：求P*Q的导函数p,q=polyder(P,Q)：求P/Q的导函数，导数分子存入p,分母存入q6.1 多项式(polynomial)Application of Matlab Language2013/9/1011多项式的估值(polyval)举例：绘制p(x)=x34x27x10在-1,3段上的曲线。x=linspace(-1,3);p=1 4-7-10;v=pol
8、yval(p,x);%choose 100 data points between-1and 3.plot(x,v);title(x3+4x2-7x-10);xlabel(x)6.1 多项式(polynomial)p(x)=x34x27x10在-1,3段上的曲线:Application of Matlab Language2013/9/10126.2 函数的数值导数导数定义为：则y=f(x)的导数可近似为：h0 f(x h)f(x)(x h)(x)limdy dxdy f(x h)f(x)dx(x h)(x)这里h0它是y的有限差分除以x的有限差分。MATLAB中没有直接提供数值导数的函数，只
9、有计算向前差分的函数diff，其调用格式为：DX=diff(X)计算向量X的向前差分DX=diff(X，n)计算向量X的n阶向前差分Application of Matlab Language2013/9/1013例题Application of Matlab Language2013/9/1014设x由0,2间均匀分布的10个点组成，求sinx的1-3阶差分。命令如下：X=linspace(0,2*pi,10);Y=sin(X);DY=diff(Y)D2Y=diff(Y,2)D3Y=diff(Y,3)DY=0.64280.34200.6428D2Y=-0.1188-0.5240-0.6840
10、-0.5240-0.11880.3420-0.3008D3Y=-0.1600-0.4608-0.4052-0.16000.16000.40520.46080.30080.05560.24520.32010.24520.0556-0.16006.2 函数的数值导数（续）例：设 f(x)x3 2 x2 x 12 6 x 5 5 x 2在-3,3区间内以0.01为步长求数值导数。并画出导函 数图像。程序如下：f=inline(sqrt(x.3+2*x.2-x+12)+(x+5).(1/6)+5*x+2);8内联函数7.5%定义式求导数x=-3:0.01:3;dx=diff(f(x,3.01)/0.0
11、1;plot(x,dx)7%根据6.56-3-2-101234.55.55Application of Matlab Language2013/9/10156.3 数值积分Application of Matlab Language2013/9/1016一元函数的数值积分常用积分指令：quad和quadl。一般说来，quadl比quad更有效。具体调用格式如下：q=quadl(fun,a,b)q=quadl(fun,a,b,tol)q=quadl(fun,a,b,tol,trace)q,fcnt=quadl(fun,a,b,.)输入量fun为被积函数的句柄。输入量a,b分别是积分的下限、和上限
12、，都必须是确定的数值;前3个输入参数是调用积分指令所必须的，其他可以缺省;输入量tol是一个标量，控制绝对误差;输入量trace为非0值时，将随积分的进程逐点画出被积分函数;输出参数fcnt返回函数的执行次数。Note：quad的调用格式与quadl相同6.3 数值积分（续）Application of Matlab Language2013/9/1017举例：求定积分 I e x dxMATLAB指令quad和quadl求积分fun=inline(exp(-x.*x),x);必须的%数组乘符号.*的采用是Isim=quad(fun,0,1),I8=quadl(fun,0,1)Isim=0.7
13、468I8=0.74682106.3 数值积分（续）举例：求解定积分 I 1Application of Matlab Language2013/9/1018用quad指令求积分ff=inline(sqrt(log(1./x),x);Isim=quad(ff,0,1)Warning:Divide by zero.In inlineeval at 13 In inline.subsref at 25 In quad at 63Isim=0.88620ln 1 dx x6.4 元素排序Application of Matlab Language2013/9/1019Matlab中对向量X排序的函数
14、是sort(X),函数返回一个 对X中的元素按升序排列的新向量。sort函数也可以对矩阵A的各列（或行）重新排序，其 调用格式为：Y,I=sort(A,dim)dim=1,按列排序；dim=2,按行排序，Y是排序后的矩阵，I记录Y中13713的元素在A中的位置。185 例：对下列矩阵做各种排序。A 4126 命令如下：A=1,-8,5;4,12,6;13,7,-13;sort(A)ans=Application of Matlab Language2013/9/1020-sort(-A,2)%对A的每行按降序排列ans=1-8-134751312651-81264137-136.5 数据插值A
15、pplication of Matlab Language2013/9/1021在工程测量和科学实验中，所得到的数据通常是离散的，要得到 这些离散点以外的其他点的数值，就需要根据已知的数据进行插值。插值函数一般由线性函数、多项式、样条函数或这些函数的分段函 数充当。一维数据插值：被插值函数有一个单变量。采用的方法有：线性方法、最近方法、三次样条和三次插值。在Matlab中实现这些插值的函数是interp1，其调用格式如下：Y1=interp1(X,Y,X1,method）函数根据X，Y的值，计算函数在X1处的值。X,Y是两个等长的已知向量，分别描述采样点和样本值；X1是一个向量或标量，描述欲插
16、值的点；Y1是一个与X1等长的插值结果。method是插值方法，允许的取值为：6.5 数据插值Application of Matlab Language2013/9/1022（1）linear：线性插值。默认的插值方式。它是把插值点靠近的两个 数据点用直线连接，然后在直线上选取对应插值点的数据。（2）nearest：最近点插值。根据已知插值点与已知数据点的远近程 度进行插值。插值点优先选择较近的数据点进行插值。（3）cubic：3次多项式插值。根据已知数据求出一个3次多项式，然 后根据该多项式进行插值。（4）spline：3次样条插值。指在每个分段内构造一个3次多项式，使 其满足插值条件外，
17、在各节点处具有光滑的条件。例：给出概率积分数据表如下，用不同的插值方法计算f(0.472)。x0.460.470.480.49f(x)0.48465550.49375420.50274980.5116683命令如下：Application of Matlab Language2013/9/1023x=0.46:0.01:0.49;f=0.4846555,0.4937542,0.5027498,0.5116683;format long interp1(x,f,0.472)ans=0.49555332000000interp1(x,f,0.472,nearest)ans=0.4937542000
18、0000interp1(x,f,0.472,spline)ans=0.49556073600000interp1(x,f,0.472,cubic)ans=0.49556111971206其中，3次样条和3次多项式的插值结果优于 最近点插值方法和线性插值方法，但插值方法 的好坏依赖于被插值函数，没有一种对所有函 数都是最好的插值方法。6.6 曲线拟合Application of Matlab Language2013/9/1024数值插值要求逼近函数在采样点与被逼近函数相等，但由于测量 误差，所获得的数据不一定准确，如果强求逼近显然不够合理。曲线拟合不要求逼近函数通过各采样点，但要尽量的接近这些
19、点，使 误差在某种意义上达到最小。曲线拟合的实现：在matlab中，用polyfit函数来求得最小二乘拟合多项式的系数，再用 polyval函数按所得的多项式计算所给出点上的函数近似值。polyfit函数的调用格式为：P,S=polyfit(X,Y,m)函数根据采样点X和采样点函数值Y，产生一个m次多项式P及其在采样 点的误差向量S。其中X、Y是两个等长的向量，P是一个长度为m+1的 向量，P的元素是多项式系数。polyval函数的功能是按多项式的系数计算x点多项式的值。例：用一个三次多项式在区间0,2内逼近函数sinx。在给定区间内，均匀的选择20个采样点，并计算采样点的函数值然后利用3次多项式逼近。1命令如下：x=linspace(0,2*pi,20);y=sin(x);p=polyfit(x,y,3)y1=polyval(p,x)plot(x,y,:o,x,y1,-*)legend(sin(x),fit)01234567-1-0.8-0.6-0.4-0.200.20.4 0.60.8 sin(x)fitApplication of Matlab Language2013/9/1025
