实验1帮助文件插值和拟合

1、实验一 插值和拟合实验一 函数的插值和拟合方法1.1 插值问题及其误差 与插值有关的MATLAB 函数(一) POLY2SYM 函数调用格式一：poly2sym (C)调用格式二：f1=poly2sym(C,'V') 或 f2=poly2sym(C, sym ('V') ),(二) POLYVAL 函数调用格式：Y = polyval(P,X)(三) POLY 函数调用格式：Y = poly (V)(四) CONV 函数调用格式：C =conv (A, B)例 1 求三个一次多项式、和的积.它们的零点分别依次为0.4,0.8,1.2.解 我们可以用两种MATLA

2、B程序求之.方法1 如输入MATLAB程序>> X1=0.4,0.8,1.2; l1=poly(X1), L1=poly2sym (l1)运行后输出结果为l1 = 1.0000 -2.4000 1.7600 -0.3840L1 = x3-12/5*x2+44/25*x-48/125方法2 如输入MATLAB程序>> P1=poly(0.4);P2=poly(0.8);P3=poly(1.2); C =conv(conv(P1, P2),P3) , L1=poly2sym (C)运行后输出的结果与方法1相同.(五) DECONV 函数调用格式：Q,R =deconv (B

3、,A)(六) roots(poly(1:n)命令调用格式：roots(poly(1:n) (七) det(a*eye(size (A) - A)命令调用格式：b=det(a*eye(size (A) - A)1.2 拉格朗日（Lagrange）插值及其MATLAB程序一线性插值及其MATLAB程序例2 已知函数在上具有二阶连续导数，且满足条件.求线性插值多项式和函数值，并估计其误差.解 输入程序>> X=1,3;Y=1,2; l01= poly(X(2)/( X(1)- X(2), l11= poly(X(1)/( X(2)- X(1), l0=poly2sym (l01),l1=

4、poly2sym (l11), P = l01* Y(1)+ l11* Y(2), L=poly2sym (P),x=1.5; Y = polyval(P,x)运行后输出基函数l0和l1及其插值多项式的系数向量P（略）、插值多项式L和插值Y为l0 = l1 = L = Y =-1/2*x+3/2 1/2*x-1/2 1/2*x+1/2 1.2500输入程序>> M=5;R1=M*abs(x-X(1)* (x-X(2)/2运行后输出误差限为 R1 = 1.8750例3 求函数e在上线性插值多项式，并估计其误差.解 输入程序>> X=0,1; Y =exp(-X) , l0

5、1= poly(X(2)/( X(1)- X(2), l11= poly(X(1)/( X(2)- X(1), l0=poly2sym (l01),l1=poly2sym (l11), P = l01* Y(1)+ l11* Y(2), L=poly2sym (P),运行后输出基函数l0和l1及其插值多项式的系数向量P和插值多项式L为l0 = l1 = P =-x+1 x -0.6321 1.0000L =-1423408956596761/2251799813685248*x+1 输入程序>> M=1;x=0:0.001:1; R1=M*max(abs(x-X(1).*(x-X(

6、2)./2运行后输出误差限为 R1 = 0.1250.二抛物线插值及其MATLAB程序例4 求将区间 0, /2 分成等份，用产生个节点，然后分别作线性插值函数和抛物线插值函数.用它们分别计算cos (/6) (取四位有效数字)，并估计其误差.解 (1)线性插值 输入程序>> X=0,pi/2; Y =cos(X) ,l01= poly(X(2)/( X(1)- X(2), l11= poly(X(1)/( X(2)- X(1), l0=poly2sym (l01),l1=poly2sym (l11), P = l01* Y(1)+ l11* Y(2), L=poly2sym (P

7、),x=pi/6; Y = polyval(P,x)运行后输出基函数l0和l1及其插值多项式的系数向量P、插值多项式和插值为l0 =-5734161139222659/9007199254740992*x+1l1 =5734161139222659/9007199254740992*xP = -0.6366 1.0000L =-5734161139222659/9007199254740992*x+1Y =0.6667输入程序>> M=1;x=pi/6; R1=M*abs(x-X(1)*(x-X(2)/2运行后输出误差限为R1 = 0.2742.（2）抛物线插值 输入程序>&

8、gt; X=0:pi/4:pi/2; Y =cos(X) ,l01= conv (poly(X(2),poly(X(3)/( X(1)- X(2)* ( X(1)- X(3), l11= conv (poly(X(1), poly(X(3)/( X(2)- X(1)* ( X(2)- X(3),l21= conv (poly(X(1), poly(X(2)/( X(3)- X(1)* ( X(3)- X(2),l0=poly2sym (l01),l1=poly2sym (l11),l2=poly2sym (l21),P = l01* Y(1)+ l11* Y(2) + l21* Y(3), L

9、=poly2sym (P),x=pi/6; Y = polyval(P,x)运行后输出基函数l01、l11和l21及其插值多项式的系数向量P、插值多项式L和插值Y为l0 =228155022448185/281474976710656*x2-2150310427208497/1125899906842624*x+1l1 =-228155022448185/140737488355328*x2+5734161139222659/2251799813685248*xl2 =228155022448185/281474976710656*x2-5734161139222659/90071992547

10、40992*xP = -0.3357 -0.1092 1.0000L=-6048313895780875/18014398509481984*x2-7870612110600739/72057594037927936*x+1Y = 0.8508输入程序>> M=1;x=pi/6; R2=M*abs(x-X(1)*(x-X(2) *(x-X(3)/6运行后输出误差限为R2 =0.0239.三 次拉格朗日（Lagrange）插值及其MATLAB程序例5 给出节点数据，作三次拉格朗日插值多项式计算，并估计其误差.解 输入程序>> X=-2,0,1,2; Y =17,1,2,1

11、7;p1=poly(X(1); p2=poly(X(2);p3=poly(X(3); p4=poly(X(4); l01= conv ( conv (p2, p3), p4)/( X(1)- X(2)* ( X(1)- X(3) * ( X(1)- X(4), l11= conv ( conv (p1, p3), p4)/( X(2)- X(1)* ( X(2)- X(3) * ( X(2)- X(4),l21= conv ( conv (p1, p2), p4)/( X(3)- X(1)* ( X(3)- X(2) * ( X(3)- X(4),l31= conv ( conv (p1, p

12、2), p3)/( X(4)- X(1)* ( X(4)- X(2) * ( X(4)- X(3),l0=poly2sym (l01),l1=poly2sym (l11),l2=poly2sym (l21), l3=poly2sym (l31),P = l01* Y(1)+ l11* Y(2) + l21* Y(3) + l31* Y(4),运行后输出基函数l0，l1，l2和l3及其插值多项式的系数向量P（略）为l0 =-1/24*x3+1/8*x2-1/12*x，l1 =1/4*x3-1/4*x2-x+1l2 =-1/3*x3+4/3*x，l3 =1/8*x3+1/8*x2-1/4*x输入程

13、序>> L=poly2sym (P),x=0.6; Y = polyval(P,x)运行后输出插值多项式和插值为L = Y =x3+4*x2-4*x+1 0.2560.输入程序>> syms M; x=0.6; R3=M*abs(x-X(1)*(x-X(2) *(x-X(3) *(x-X(4)/24运行后输出误差限为R3 =91/2500*M即 R3 ， .四拉格朗日多项式和基函数的MATLAB程序求拉格朗日插值多项式和基函数的MATLAB主程序function C, L,L1,l=lagran1(X,Y)m=length(X); L=ones(m,m);for k=1

14、: m V=1; for i=1:m if k=i V=conv(V,poly(X(i)/(X(k)-X(i); endendL1(k,:)=V; l(k,:)=poly2sym (V)， endC=Y*L1;L=Y*l例6 给出节点数据， ，作五次拉格朗日插值多项式和基函数，并写出估计其误差的公式.解 在MATLAB工作窗口输入程序>> X=-2.15 -1.00 0.01 1.02 2.03 3.25;Y=17.03 7.24 1.05 2.03 17.06 23.05;C, L ,L1,l= lagran1(X,Y)运行后输出五次拉格朗日插值多项式L及其系数向量C，基函数l及

15、其系数矩阵L1如下C = -0.2169 0.0648 2.1076 3.3960 -4.5745 1.0954L =1.0954-4.5745*x+3.3960*x2+2.1076*x3+0.0648*x4-0.2169*x5L1 = -0.0056 0.0299 -0.0323 -0.0292 0.0382 -0.0004 0.0331 -0.1377 -0.0503 0.6305 -0.4852 0.0048 -0.0693 0.2184 0.3961 -1.2116 -0.3166 1.0033 0.0687 -0.1469 -0.5398 0.6528 0.9673 -0.0097

16、-0.0317 0.0358 0.2530 -0.0426 -0.2257 0.0023 0.0049 0.0004 -0.0266 0.0001 0.0220 -0.0002l = -0.0056*x5+0.0299*x4-0.0323*x3-0.0292*x2+0.0382*x-0.0004 0.0331*x5-0.1377*x4-0.0503*x3+0.6305*x2-0.4852*x+0.0048 -0.0693*x5+0.2184*x4+0.3961*x3-1.2116*x2-0.3166*x+1.0033 0.0687*x5-0.1469*x4-0.5398*x3+0.6528*x

17、2+0.9673*x-0.0097 -0.0317*x5+0.0358*x4+0.2530*x3-0.0426*x2-0.2257*x+0.0023 0.0049*x5+0.0004 *x4-0.0266*x3+0.0001*x2+0.0220*x-0.0002估计其误差的公式为，.五拉格朗日插值及其误差估计的MATLAB程序拉格朗日插值及其误差估计的MATLAB主程序function y,R=lagrange(X,Y,x,M)%X为采样点，Y为采样点的函数值，x为插值点，M为函数n阶导数的最大值n=length(X); m=length(x);for i=1:m z=x(i);s=0.0;

18、for k=1:n p=1.0; q1=1.0; c1=1.0;for j=1:n if j=kp=p*(z-X(j)/(X(k)-X(j);%基函数 end q1=abs(q1*(z-X(j);c1=c1*j;%插值余项 end s=p*Y(k)+s; end y(i)=s;%最终的结果endR=M*q1/c1;%最终的余项值例 7 已知，用拉格朗日插值及其误差估计的MATLAB主程序求的近似值，并估计其误差.解 在MATLAB工作窗口输入程序>> x=2*pi/9; M=1; X=pi/6 ,pi/4, pi/3;Y=0.5,0.7071,0.8660; y,R=lagrang

19、e(X,Y,x,M)运行后输出插值y及其误差限R为y = R =0.6434 8.8610e-004.1.2 曲线拟合、误差及其MATLAB程序例8 已知函数和一组数据列入表1中，比较最大误差，平均误差，均方根误差和误差平方和. 表1 例8的一组数据xi-2.5 -1.7 -1.1 -0.8 0 0.1 0.5 3.6yi-43.50 5.69 11.34 14.16 0 1.02 -6.37 185.84解 由给定的函数和数据，在MATLAB工作窗口输入>> x=-2.5,-1.7,-1.1,-0.8,0,0.1,0.5,3.6; n=length(x);y=-43.50 5.6

20、9 11.34 14.16 0 1.02 -6.37 185.84;f=5.*x.3-14.*x+7.*(sin(2*pi*x).2; fy=abs(f-y);fy2=fy.2; x',y',f',fy',fy2', Ew=max(fy),E1=sum(fy)/n, E2=sqrt(sum(fy2)/n), E=sum(fy2)运行后屏幕显示如下 x y f fy fy2-2.5000 -43.5000 -43.1250 0.3750 0.1406-1.7000 5.6900 5.5666 0.1234 0.0152-1.1000 11.3400 11.

21、1634 0.1766 0.0312-0.8000 14.1600 14.9716 0.8116 0.6586 0 0 0 0 0 0.1000 1.0200 1.0234 0.0034 0.0000 0.5000 -6.3700 -6.3750 0.0050 0.0000 3.6000 185.8400 185.2984 0.5416 0.2933Ew = E1 = E2 = E =0.8116 0.2546 0.3773 1.1390曲线拟合的线性最小二乘法及其MATLAB程序例9 给出一组数据点列入表2中，试用线性最小二乘法求拟合曲线，并用估计其误差，作出拟合曲线.表2 例8 的一组数据

22、xi-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6yi-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04解 （1）在MATLAB工作窗口输入程序>> x=-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6; y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04;plot(x,y,'r*'),legend('实验数据(xi,yi)')xlabel('x'),

23、 ylabel('y'),title('例8的数据点(xi,yi)的散点图')运行后屏幕显示数据的散点图（略）. （2）编写下列MATLAB程序计算在处的函数值，假设拟合函数为三次多项式，输入程序>> syms a1 a2 a3 a4x=-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6; fi=a1.*x.3+ a2.*x.2+ a3.*x+ a4运行后屏幕显示关于a1,a2, a3和a4的线性方程组，即为超越方程组fi = -125/8*a1+25/4*a2-5/2*a3+a4, -4913/1000*a1+289/100

24、*a2-17/10*a3+a4, -1331/1000*a1+121/100*a2-11/10*a3+a4, -64/125*a1+16/25*a2-4/5*a3+a4, a4, 1/1000*a1+1/100*a2+1/10*a3+a4, 27/8*a1+9/4*a2+3/2*a3+a4, 19683/1000*a1+729/100*a2+27/10*a3+a4, 5832/125*a1+324/25*a2+18/5*a3+a4编写构造误差平方和的MATLAB程序>> y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68

25、.04;fi = -125/8*a1+25/4*a2-5/2*a3+a4, -4913/1000*a1+289/100*a2-17/10*a3+a4, -1331/1000*a1+121/100*a2-11/10*a3+a4, -64/125*a1+16/25*a2-4/5*a3+a4, a4, 1/1000*a1+1/100*a2+1/10*a3+a4, 27/8*a1+9/4*a2+3/2*a3+a4, 19683/1000*a1+729/100*a2+27/10*a3+a4, 5832/125*a1+324/25*a2+18/5*a3+a4fy=fi-y; fy2=fy.2; J=sum

26、(fy.2)%误差平方和运行后屏幕显示误差平方和如下J= (-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)2+(-4913/1000*a1+289/100*a2-17/10*a3+a4+171/2)2+(-1331/1000*a1+121/100*a2-11/10*a3+a4+723/20)2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)2+(a4+91/10)2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)2+(19683/1000

27、*a1+729/100*a2+27/10*a3+a4-13/2)2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)2为求使达到最小，只需利用极值的必要条件 ，得到关于的线性方程组，这可以由下面的MATLAB程序完成，即输入程序 >> syms a1 a2 a3 a4J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)2+(-4913/1000*a1+289/100*a2-17/10*a3+a4.+171/2)2+(-1331/1000*a1+121/100*a2-11/10*a3+a4+723/20)2+(-64/12

28、5*a1+16/25*a2-4/5*a3+a4+663/25)2+(a4+91/10)2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/2)2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)2; Ja1=diff(J,a1); Ja2=diff(J,a2); Ja3=diff(J,a3); Ja4=diff(J,a4);Ja11=simple(Ja1), Ja21=sim

29、ple(Ja2), Ja31=simple(Ja3), Ja41=simple(Ja4),运行后屏幕显示J分别对a1, a2 ,a3 ,a4的偏导数如下Ja11=56918107/10000*a1+32097579/25000*a2+1377283/2500*a3+23667/250*a4-8442429/625Ja21 = 32097579/25000*a1+1377283/2500*a2+23667/250*a3+67*a4+767319/625Ja31 =1377283/2500*a1+23667/250*a2+67*a3+18/5*a4-232638/125Ja41 =23667/25

30、0*a1+67*a2+18/5*a3+18*a4+14859/25解线性方程组Ja11 =0，Ja21 =0，Ja31 =0，Ja41 =0，输入下列程序>>A=56918107/10000, 32097579/25000, 1377283/2500, 23667/250; 32097579/25000, 1377283/2500, 23667/250, 67; 1377283/2500, 23667/250, 67, 18/5; 23667/250, 67, 18/5, 18;B=8442429/625, -767319/625, 232638/125, -14859/25; C

31、=B/A, f=poly2sym(C)运行后屏幕显示拟合函数f及其系数C如下C = 5.0911 -14.1905 6.4102 -8.2574f=716503695845759/140737488355328*x3-7988544102557579/562949953421312*x2+1804307491277693/281474976710656*x-4648521160813215/562949953421312 故所求的拟合曲线为.（3）编写下面的MATLAB程序估计其误差，并作出拟合曲线和数据的图形.输入程序>> xi=-2.5 -1.7 -1.1 -0.8 0 0.1

32、 1.5 2.7 3.6; y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04;n=length(xi); f=5.0911.*xi.3-14.1905.*xi.2+6.4102.*xi -8.2574;x=-2.5:0.01: 3.6; F=5.0911.*x.3-14.1905.*x.2+6.4102.*x -8.2574;fy=abs(f-y); fy2=fy.2; Ew=max(fy), E1=sum(fy)/n, E2=sqrt(sum(fy2)/n)plot(xi,y,'r*'), hold on,

33、 plot(x,F,'b-'), hold offlegend('数据点(xi,yi)','拟合曲线y=f(x)'), xlabel('x'), ylabel('y'),title('例9的数据点(xi,yi)和拟合曲线y=f(x)的图形')运行后屏幕显示数据与拟合函数f的最大误差Ew，平均误差E1和均方根误差E2及其数据点和拟合曲线y=f(x)的图形（略）.Ew = E1 = E2 =3.105 4 0.903 4 1.240 91.3 多项式拟合及其MATLAB程序例7.4.1 给出一组数据点列入表73中，试用线性最小二乘法求拟合曲线，并用（7.2），（7.3）和（7.4）式估计其误差，作出拟合曲线.表73 例7