MATLAB插 值 法

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