《matlab应用(多项式插值)》由会员分享，可在线阅读，更多相关《matlab应用(多项式插值)（54页珍藏版）》请在金锄头文库上搜索。
1、二 多项式与插值 多项式插值的主要目的是用一个多项式 拟合离散点上的函数值，使得可以用该 多项式估计数据点之间的函数值。 可导出数值积分方法，有限差分近似 关注插值多项式的表达式、精度、选点 效果。2.1 多项式MATLAB命令 一个多项式的幂级数形式可表示为： 也可表为嵌套形式 或因子形式N阶多项式n个根，其中包含重根和复根。若多 项式所有系数均为实数，则全部复根都将以共轭对 的形式出现 幂系数：在MATLAB里，多项式用行向量表示，其 元素为多项式的系数，并从左至右按降幂排列。例：被表示为 p=2 1 4 5 poly2sym(p)ans =2*x3+x2+4*x+5 Roots: 多项式
2、的零点可用命令roots求的。例： r=roots(p) 得到r =0.2500 + 1.5612i0.2500 - 1.5612i-1.0000 所有零点由一个列向量给出。 Poly: 由零点可得原始多项式的各系数，但可能相差 一个常数倍。例： poly(r) ans =1.0000 0.5000 2.0000 2.5000注意：若存在重根，这种转换可能会降低精度。例： r=roots(1 -6 15 -20 15 -6 1) r =1.0042 + 0.0025i1.0042 - 0.0025i1.0000 + 0.0049i1.0000 - 0.0049i0.9958 + 0.0024i
3、0.9958 - 0.0024i 舍入误差的影响，与计算精度有关。 polyval: 可用命令polyval计算多项式的值 。例： 计算y(2.5) c=3,-7,2,1,1； xi=2.5; yi=polyval(c,xi)yi =23.8125 如果xi是含有多个横坐标值的数组，则yi也 为与xi长度相同的向量。 c=3,-7,2,1,1; xi=2.5,3; yi=polyval(c,xi) yi =23.8125 76.0000 polyfit:给定n+1个点将可以唯一确定一个n阶多项式。利 用命令polyfit可容易确定多项式的系数。例： x=1.1,2.3,3.9,5.1; y=3
4、.887,4.276,4.651,2.117; a=polyfit(x,y,length(x)-1) a =-0.2015 1.4385 -2.7477 5.4370 poly2sym(a)ans =-403/2000*x3+2877/2000*x2-27477/10000*x+5437/1000 多项式为 Polyfit的第三个参数是多项式的阶数。多项式积分：功能：求多项式积分调用格式：py=poly_itg(p)p:被积多项式的系数py：求积后多项式的系数poly_itg.mfunction py=poly_itg(p)n=length(p);py=p.*n:-1:1.(-1),0 不包括
5、最后一项积分常数多项式微分： Polyder: 求多项式一阶导数的系数。调用格式为: b=polyder(c )c为多项式y的系数，b是微分后的系数，其值为：两个多项式的和与差：命令poly_add:求两个多项式的和,其调用格式为:c= poly_add(a,b)多项式a减去b，可表示为:c= poly_add(a,-b)功能：两个多项式相加调用格式：b=poly_add(p1,p2)b:求和后的系数数组 poly_add.m function p3=poly_add(p1,p2) n1=length(p1); n2=length(p2); if n1=n2 p3=p1+p2;end if n
6、1n2 p3=p1+zeros(1,n1-n2),p2;end if n1 a=2,-5,6,-1,9; b=3,-90,-18; c=conv(a,b) c =6 -195 432 -453 9 -792 -162 q,r=deconv(c,b) q =2 -5 6 -1 9 r =0 0 0 0 0 0 0 poly2sym(c)ans =6*x6-195*x5+432*x4-453*x3+9*x2-792*x-1622.2 Lagrange插值 方法介绍对给定的n个插值点 及对应的函 数值 ，利用n次Lagrange插值多项式 ，则对插值区间内任意x的函数值y可通过下 式求的： MATL
7、AB实现function y=lagrange(x0,y0,x) ii=1:length(x0); y=zeros(size(x); for i=iiij=find(ii=i); y1=1;for j=1:length(ij), y1=y1.*(x-x0(ij(j); endy=y+y1*y0(i)/prod(x0(i)-x0(ij); end 算例：给出f(x)=ln(x)的数值表，用Lagrange计算 ln(0.54)的近似值。 x=0.4:0.1:0.8； y=-0.916291,-0.693147,-0.510826,-0.356675,- 0.223144; lagrange(x,
8、y,0.54) ans =-0.6161 （精确解-0.616143)2.3 Hermite插值 方法介绍不少实际问题不但要求在节点上函数值相等，而且 要求导数值也相等，甚至要求高阶导数值也相等，满足 这一要求的插值多项式就是Hermite插值多项式。下面 只讨论函数值与一阶导数值个数相等且已知的情况。已知n个插值点 及对应的函数值和一阶导数值 。则对插值区间 内任意x的函数值y的Hermite插值公式： MATLAB实现 % hermite.m function y=hermite(x0,y0,y1,x) n=length(x0); m=length(x); for k=1:m yy=0.0
9、;for i=1:n h=1.0; a=0.0;for j=1:nif j=ih=h*(x(k)-x0(j)/(x0(i)-x0(j)2;a=1/(x0(i)-x0(j)+a;endendyy=yy+h*(x0(i)-x(k)*(2*a*y0(i)-y1(i)+y0(i);endy(k)=yy; end 算例:对给定数据,试构造Hermite多项式求出 sin0.34的近似值。 x0=0.3,0.32,0.35; y0=0.29552,0.31457,0.34290; y1=0.95534,0.94924,0.93937; y=hermite(x0,y0,y1,0.34) y =0.3335
10、sin(0.34) 与精确值比较 ans =0.3335 x=0.3:0.005:0.35;y=hermite(x0,y0,y1,x); plot(x,y) y2=sin(x); hold on plot(x,y2,-r)2.4 Runge现象和分段插值 问题的提出：根据区间a,b上给出的节点做 插值多项式p(x)的近似值，一般总认为p(x)的 次数越高则逼近f(x)的精度就越好，但事实并 非如此。 反例：在区间-5,5上的各阶导数存在，但在此 区间上取n个节点所构成的Lagrange插值多项 式在全区间内并非都收敛。 取n=10,用Lagrange插值法进行插值计算。 x=-5:1:5; y
11、=1./(1+x.2); x0=-5:0.1:5; y0=lagrange(x,y,x0); y1=1./(1+x0.2);绘制图形 plot(x0,y0,-r)插值曲线 hold on plot(x0,y1,-b)原曲线 为解决Rung问题，引入分段插值。 算法分析：所谓分段插值就是通过插值点用折 线或低次曲线连接起来逼近原曲线。 MATLAB实现 可调用内部函数。 命令1 interp1 功能 ： 一维数据插值（表格查找）。该命令对数据点之 间计算内插值。它找出一元函数f(x)在中间点的数值。其 中函数f(x)由所给数据决定。 格式1 yi = interp1(x,Y,xi) %返回插值向
12、量yi，每一元素对应于参量xi，同时由向 量x与Y的内插值决定。参量x指定数据Y的点。若Y为一矩阵 ，则按Y的每列计算。 算例 对于t，beta 、alpha分别有两组数据与之对应，用分段线 性插值法计算当t=321, 440, 571时beta 、alpha的值。2.5 分段插值 temp=300,400,500,600; beta=1000*3.33,2.50,2.00,1.67; alpha=10000*0.2128,0.3605,0.5324,0.7190; ti=321,400,571; propty=interp1(temp,beta,alpha,ti); propty=inter
13、p1(temp,beta,alpha,ti ,linear); ti,propty ans =1.0e+003 *0.3210 3.1557 2.43820.4000 2.5000 3.60500.5710 1.7657 6.6489 格式2 yi = interp1(Y,xi) %假定x=1:N，其中N为向量Y的长度，或者为矩阵 Y的行数。 格式3 yi = interp1(x,Y,xi,method) %用指定的算法计算 插值： nearest：最近邻点插值，直接完成计算；linear：线性插值（缺省方式），直接完成计算 ；spline：三次样条函数插值。cubic： 分段三次Hermit
14、e插值。其它，如v5cubic 。对于超出x范围的xi的分量，使用方法nearest、 linear、v5cubic的插值算法，相应地将返回NaN。 对其他的方法，interp1将对超出的分量执行外插值算 法。 yi = interp1(x,Y,xi,method,extrap) yi = interp1(x,Y,xi,method,extrapval) %确定超出x范围 的xi中的分量的外插值extrapval，其值通常取NaN或0 。 算例 year=1900:10:2010; product=75.995,91.972,105.711,123.203,131.669,. 150.697,179.323,203.212,226.505,249.633,256.344,267.893; p1995 = interp1(year,product,1995) p1995 =252.9885 x = 1900:1:2010; y = interp1(year, product,x,cubic); plot(year, produc
