<PRE class=m name="code">%1.n個節點分段Lagrange插值多項式；   
%2.使用格式y=lagrange(x0,y0,x,k)；   
%3.輸入項x0為n維插值節點向量，y0為n維被插函數值向量；   
%4.x為m維插值點向量，k為分段插值多項式次數，不超過3，缺省為k=1；   
%5.輸出y為插值點x處的函數值；   
%6.本程序于2002.4.21.編寫?    
function y=lagrange(x0,y0,x,k)   
if nargin<4,k=1;end   
if k>3,error('分段次數過高，容易產生Runge現象，請重新選擇次數k.'),end   
n=length(x0);   
m=length(x);   
nn=1;   
for i=1:m   
   u=x(i);   
   switch k   
   %---------------------------針對不同的k判斷插值區間      
   case 1   
      if u<=x0(2)   
         t=1;     
      elseif u>x0(n-1)   
         t=n-1;    
      else   
         for j=nn+1:n-2   
            if u>x0(j)&u<=x0(j+1)   
               t=j;nn=t-1;break   
            end   
         end   
      end   
   %---------------------------      
   case 2   
      if u<=(x0(2)+x0(3))/2   
         t=1;   
      elseif u>(x0(n-2)+x0(n-1))/2   
         t=n-2;   
      else   
         for j=nn+1:n-3   
            if u>(x0(j)+x0(j+1))/2 & u<=(x0(j+1)+x0(j+2))/2   
               t=j;nn=t-1;break   
            end   
         end   
      end   
   %---------------------------      
   case 3   
      if u<=x0(3)   
         t=1;   
      elseif u>x0(n-2)   
         t=n-3;   
      else   
         for j=nn+2:n-3   
            if u>x0(j) & u<=x0(j+1)   
               t=j-1;nn=t;break   
            end   
         end   
      end   
   end   
   %--------------------------------主程序   
  s=0;   
  for j=t:t+k   
      L=1;   
       for p=t:t+k   
          if p~=j   
            L=L*(u-x0(p))/(x0(j)-x0(p));   
          end   
       end   
      s=s+L*y0(j);   
    end   
   y(i)=s;   
end   
</PRE>  
<SCRIPT src="/inc/gg_read2.js"></SCRIPT>