poly_d2l.m

这是国外用的研究分岔的完整的M程序

function d2p=poly_d2l(m,c);

% function d2p=poly_d2l(m,c);
% INPUT:
%       m lagrange degree
%       c evaluation point in [0,1]
% OUTPUT:
%       d2p values of second derivative of lagrange polynomials 
%           through 0:1/m:1 at c

% (c) DDE-BIFTOOL v. 2.00, 30/11/2001 

if m==1
  d2p(1)=0;
  d2p(2)=0;
elseif m==2
  d2p(1)=4;
  d2p(2)=-8;
  d2p(3)=4;
elseif m==3
  d2p(1)=-13.5*(c-0.8591167563965418+c-0.4742165769367915);
  d2p(2)=40.5*(c-0.8495279234516212+c-0.2615831876594899);
  d2p(3)=-40.5*(c-0.7384168123405100+c-0.1504720765483788);
  d2p(4)=13.5*(c-0.5257834230632086+c-0.1408832436034581);
elseif m==4
  d2p(1)=42.666666666666667*((c-0.625)*(c-0.9045084971874737)+ ...
               (c-0.3454915028125263)*(c-0.9045084971874737)+ ...
               (c-0.3454915028125263)*(c-0.625));
  d2p(2)=-170.666666666666667*((c-0.61727336112119537025)* ... 
                               (c-0.90181300471826281827)+ ...
               (c-0.16841363416054181147)*(c-0.90181300471826281827)+ ...
               (c-0.16841363416054181147)*(c-0.61727336112119537025));
  d2p(3)=256*((c-0.5)*(c-0.8952847075210475)+ ...
              (c-0.1047152924789526)*(c-0.8952847075210475)+ ...
              (c-0.1047152924789526)*(c-0.5));
  d2p(4)=-170.666666666666667*((c-0.38272663887880462974)* ...
                               (c-0.83158636583945818852)+ ...
              (c-0.09818699528173718172)*(c-0.83158636583945818852)+ ...
              (c-0.09818699528173718172)*(c-0.38272663887880462974));
  d2p(5)=42.666666666666667*((c-0.375)*(c-0.6545084971874737)+ ...
              (c-0.09549150281252627)*(c-0.6545084971874737)+ ...
              (c-0.09549150281252627)*(c-0.375));
else
  t=0:1/m:1;
  for j=1:m+1
    d2p(j)=0;
    for k=1:m+1
      if k~=j
        ss=0;
        for s=1:m+1
          if s~=j & s~=k
            f=1;
            for l=1:m+1
              if l~=k & l~=j & l~=s
                f=f*(c-t(l))/(t(j)-t(l));  
              end;
            end;
            ss=ss+f/(t(j)-t(s));
          end;
        end;
        d2p(j)=d2p(j)+ss/(t(j)-t(k));
      end;
    end;
  end;
end;

return;