MATLAB插 值 法
作者:凯鲁嘎吉 - 博客园
http://www.cnblogs.com/kailugaji/
一、实验目的
二、实验原理
三、实验程序
四、实验内容
五、解答
1. 程序
(1)Lagrange插值多项式
function [C, L,L1,l]=lagran1(X,Y) %输出C为插值多项式的系数,L为插值多项式,L1为l的系数,l为小l %输入数据表X=[];Y=[];行向量 m=length(X); L=ones(m,m); for k=1: m V=1; for i=1:m if k~=i V=conv(V,poly(X(i)))/(X(k)-X(i)); end end L1(k,:)=V; l(k,:)=poly2sym (V); end C=Y*L1;L=Y*l;
(2)Newton插值多项式
function [A,C,L,wcgs,Cw]= newploy(X,Y) n=length(X); A=zeros(n,n); A(:,1)=Y'; q=1.0; c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1)); end b=poly(X(j-1));q1=conv(q,b); c1=c1*j; q=q1; end C=A(n,n); b=poly(X(n)); q1=conv(q1,b); for k=(n-1):-1:1 C=conv(C,poly(X(k))); d=length(C); C(d)=C(d)+A(k,k); end L(k,:)=poly2sym(C); Q=poly2sym(q1); syms M wcgs=M*Q/c1; Cw=q1/c1;
2. 运算结果
(1) >> X=[0:0.4:2]; >> Y=X.^4; >> [C, L,L1,l]=lagran1(X,Y) C = 0.0000 1.0000 0 -0.0000 0 0 L = x^4 L1 = -0.8138 4.8828 -11.0677 11.7188 -5.7083 1.0000 4.0690 -22.7865 46.2240 -40.1042 12.5000 0 -8.1380 42.3177 -76.8229 55.7292 -12.5000 0 8.1380 -39.0625 63.8021 -40.6250 8.3333 0 -4.0690 17.9036 -26.6927 15.8854 -3.1250 0 0.8138 -3.2552 4.5573 -2.6042 0.5000 0 l = - (625*x^5)/768 + (625*x^4)/128 - (2125*x^3)/192 + (375*x^2)/32 - (137*x)/24 + 1 (3125*x^5)/768 - (4375*x^4)/192 + (8875*x^3)/192 - (1925*x^2)/48 + (25*x)/2 - (3125*x^5)/384 + (8125*x^4)/192 - (7375*x^3)/96 + (2675*x^2)/48 - (25*x)/2 (3125*x^5)/384 - (625*x^4)/16 + (6125*x^3)/96 - (325*x^2)/8 + (25*x)/3 - (3125*x^5)/768 + (6875*x^4)/384 - (5125*x^3)/192 + (1525*x^2)/96 - (25*x)/8 (625*x^5)/768 - (625*x^4)/192 + (875*x^3)/192 - (125*x^2)/48 + x/2 (2) >> X=[0:0.4:2]; >> Y=X.^4; >> [A,C,L,wcgs,Cw]= newploy(X,Y) A = 0 0 0 0 0 0 0.0256 0.0640 0 0 0 0 0.4096 0.9600 1.1200 0 0 0 2.0736 4.1600 4.0000 2.4000 0 0 6.5536 11.2000 8.8000 4.0000 1.0000 0 16.0000 23.6160 15.5200 5.6000 1.0000 0.0000 C = 0.0000 1.0000 0.0000 -0.0000 0.0000 0 L = (57*x^5)/18014398509481984 + x^4 + (209*x^3)/9007199254740992 - (525*x^2)/36028797018963968 + (213*x)/72057594037927936 wcgs = -(M*(- x^6 + 6*x^5 - (68*x^4)/5 + (72*x^3)/5 - (4384*x^2)/625 + (768*x)/625))/720 Cw = 0.0014 -0.0083 0.0189 -0.0200 0.0097 -0.0017 0
3. 拓展
function [y,R]=lagran2(X,Y,x,M) %输入X=[];Y=[];行向量,x预测点,可以一个,也可以为矩阵x=[];M为x的个数, n=length(X); m=length(x); for i=1:m z=x(i);s=0.0; for k=1:n p=1.0; q1=1.0; c1=1.0; for j=1:n if j~=k p=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; end R=M*q1/c1;
在MATLAB工作窗口输入程序
>> x=2*pi/9; M=1; X=[pi/6 ,pi/4, pi/3];
Y=[0.5,0.7071,0.8660]; [y,R]=lagran2(X,Y,x,M)
运行后输出插值y及其误差限R为
y =
0.6434
R =
8.8610e-04