MatLab(3)

拉格朗日插值

function[y] = Lagrange(x)
a = [25 40 50 60];
b = [95 75 63 54];
n = 4;
sum = 0;
for k = 1:n
    pro = 1;
    for i = 1:n
        if i ~= k
            pro = pro * (x - a(i))/(a(k) - a(i));
        end
    end
    sum = sum + b(k) * pro;
end
y = sum;
end

 

牛顿插值法

function[y] = Newton(x)
n = 10;
T = zeros(n,n);
a = zeros(1,n);
b = [67.052 68.008 69.803 72.024 73.400 72.063 74.669 74.487 74.065 76.777];
a(1) = 1994;
for i = 2:n
    a(i) = a(i - 1) + 1;
end
for i = 1:n
    T(i,1) = b(i);
end
for i = 2:n
    for j = 2:i
        T(i,j) = (T(i-1,j-1) - T(i,j-1))/(a(i - j + 1) - a(i));
    end
end
sum = b(1);
pro = ones(n - 1);
pro(1) = x - a(1);
for i = 2:n-1
    pro(i) = pro(i - 1) * (x - a(i));
end
for i = 2:n
    sum = sum + T(i,i) * pro(i - 1);
end
y = sum;
m = [a x];
n = [b sum];
end

三次插条

x=[0,1,2,3];
y=[3,5,4,1];
y0=1;
yn=1;

M=spline_2(x,y,y0,yn);
disp(M);
M=[y0;M;yn];

x1=0:0.1:3;
s=zeros(size(x1,2),1);
a=[];
b=[];
for j=1:size(x1,2)
    for i=1:(size(x,2)-1)       
        if x(i)<=x1(j)&&x1(j)<=x(i+1)
             p1=M(i)*(x(i+1)-x1(j))^3/(6*(x(i+1)-x(i)));
             p2=M(i+1)*(x1(j)-x(i))^3/(6*(x(i+1)-x(i)));
             p3=(y(i)-M(i)*(x(i+1)-x(i))^2/6)*(x(i+1)-x1(j))/(x(i+1)-x(i));
             p4=(y(i+1)-M(i+1)*(x(i+1)-x(i))^2/6)*(x1(j)-x(i))/(x(i+1)-x(i));
            s(j,1)=p1+p2+p3+p4;
            a=[a x1(j)];
            b=[b s(j)];
        end
    end
end
hold on
grid on
plot(a,b);

function M=spline_2(X,Y,y0,yn)
n=length(X); 
A=zeros(n,n);
A(:,1)=Y';
D=zeros(n-2,n-2);
d=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

%求h_i
for i=1:n-1
    h(i)=X(i+1)-X(i);
end

%求mu_i,lamda_i
for i=2:n-2
    mu(i)=h(i)/(h(i)+h(i+1));
    lamda(i-1)=h(i)/(h(i-1)+h(i));
    D(i-1,i)=lamda(i-1);
    D(i,i-1)=mu(i);
end
 
%求d_i
for i=1:n-2
     D(i,i)=2;
     if (i==1)
         d(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))
         d(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
         d(i,1)=(6/(h(i+1)+h(i)))*(A(i+2,2)-A(i+1,2));
     end             
end
M=D\d;
end