main.m

clear;
clc;

x=[25 40 50 60];
y=[95 75 63 54];
xh=70;
lagrange(x,y,xh)

lagrange.m

function yh=lagrange(x,y,xh)

n = length(x); 
m = length(xh); 
p = zeros(n,m);

for k = 1:n
   t = ones(n,m);
   for j = 1:n
       if j~=k
           if abs(x(k) - x(j))<eps
              error('% 输入的插值节点必须互异！');
           end
           t(j,:) = (xh - x(j))/(x(k) - x(j));
       end
   p(k,:) = prod(t);
   end
   yh = y*p;
end

main.m

clear;
clc;

x0=[1994 1995 1996 1997 1998 1999 2000 2001 2002 2003];
y0=[67.052 68.008 69.803 72.024 73.400 72.063 74.669 74.487 74.065 76.777];
xi=2010;

q=Newton(xi,x0,y0)

Newton.m

function p= Newton(x,xi,yi)
n=length(xi);
f=zeros(n,n);

% 对差商表第一列赋值
for k=1:n      
    f(k)=yi(k);
end

% 求差商表
for i=2:n       % 差商表从0阶开始；但是矩阵是从1维开始存储！！！！！！
    for k=i:n
        f(k,i)=(f(k,i-1)-f(k-1,i-1))/(xi(k)-xi(k+1-i));  
    end
end
disp('差商表如下：');
disp(f);

%求插值多项式
p=0;          
for k=2:n
    t=1;
    for j=1:k-1
        t=t*(x-xi(j));
        %disp(t)
    end
    p=f(k,k)*t+p;
    %disp(p)
end
p=f(1,1)+p;

end

 

​

clear;
clc;
%第二类边界条件 M0=M3=1 在0处的二阶倒=1处的二阶导=1
x = [0 1 2 3];     %设置4个控制点
y = [3 5 4 1];
y0=0;              %  S''(x0)=f''(x0)=y0   
yn=0;              %  S''(xn)=f''(xn)=yn
x0=0:0.01:3;
s=threesimple2(x,y,x0,y0,yn)
plot(x0,s)        %绘制第二边界条件插值函数图像
hold on
grid on
plot(x,y,'o')
%axis([0.2 0.55 0.4 0.75])
xlabel('自变量 X'), ylabel('因变量 Y')
title('插值点与三次样条函数') 
legend('三次样条插值点坐标','插值点')

%方法一：使用spline函数    第二种边界条件，二阶导导数确定的
cs2 = spline(x,[1,y,1]);  % 1为二阶导的值
xx = linspace(x(1),x(end),100);
yy=ppval(cs2,xx);
figure(2)
plot(x,y,'bo',xx,yy,'r-');

%方法二：此次编写的第二边界matlab代码其中包含了自然边界

function [D,h,A,g,M]=three2(X,Y,y0,yn)
%        第二边界条件的三次样条函数(包含自然边界条件)
%        y0,yn表示的是S''(x0)=f''(x0)=y0，S''(xn)=f''(xn)=yn,自然边界即条件值为0 
%        此函数为M值求值函数
%        D,h,A,g,M输出量分别为系数矩阵D，插值宽度h，差商表A，g值,M值 
         n=length(X); 
         A=zeros(n,n);A(:,1)=Y';D=zeros(n-2,n-2);g=zeros(n-2,1);
         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
         end
         
         for i=1:n-1
             h(i)=X(i+1)-X(i);
         end        
         for i=1:n-2
             D(i,i)=2;
             if (i==1)
                 g(i,1)=(6/(h(i+1)+h(i)))*(A(i+2,2)-A(i+1,2))-h(i)/(h(i)+h(i+1))*y0;
             elseif (i==(n-2))
                 g(i,1)=(6/(h(i+1)+h(i)))*(A(i+2,2)-A(i+1,2))-(1-h(i)/(h(i)+h(i+1)))*yn;
             else
                 g(i,1)=(6/(h(i+1)+h(i)))*(A(i+2,2)-A(i+1,2));
             end             
         end
         for i=2:n-2
             u(i)=h(i)/(h(i)+h(i+1));
             n(i-1)=h(i)/(h(i-1)+h(i));
             D(i-1,i)=n(i-1);
             D(i,i-1)=u(i);             
         end
         M=D\g;
         M=[y0;M;yn];         
end

function s=threesimple2(X,Y,x,y0,yn)
%        第二边界条件函数 
%        s函数表示三次样条插值函数插值点对应的函数值
%        根据三次样条参数函数求出的D,h,A,g,M
%        x表示求解插值点函数点，X为已知插值点        
         [D,h,A,g,M]=three2(X,Y,y0,yn)
         n=length(X); m=length(x);    
         for t=1:m
            for i=1:n-1
               if (x(t)<=X(i+1))&&(x(t)>=X(i))
                  p1=M(i,1)*(X(i+1)-x(t))^3/(6*h(i));
                  p2=M(i+1,1)*(x(t)-X(i))^3/(6*h(i));
                  p3=(A(i,1)-M(i,1)/6*(h(i))^2)*(X(i+1)-x(t))/h(i);
                  p4=(A(i+1,1)-M(i+1,1)/6*(h(i))^2)*(x(t)-X(i))/h(i);
                  s(t)=p1+p2+p3+p4; 
                  break;
               else
                   s(t)=0; 
               end
            end
         end
end