数值计算插值

hermit三次插值

function  s_=her_three(x,y,m,x_value)
%% 输入值分配，x_input,y_input均为数组,y_0,y_n为x_0,x_n分别对应的一阶导数值
x_input = x;
y_input = y;
%% 
number = length(x_input);
delta_h = zeros(1,number-1); %给delta_h分配数组大小1×(number-1)，并全部初始化为0
delta_f = zeros(1,number-1); %给delta_f分配数组大小1×(number-1)，并全部初始化为0
 %给e分配数组大小1×(number-2)，并全部初始化为0
% 计算delta_h、delta_f的值
for i = 1:(number-1)
    delta_h(i) = x_input(i+1) - x_input(i);
    delta_f(i) = (y_input(i+1) - y_input(i))/ delta_h(i);
end


for i =1:number-1
    % 获取相邻两点间的插值函数
    x_ = linspace(x_input(1,i),x_input(1,i+1));%这里的x_相当于sym x，但是取了很多点
%     disp('x_');
%     disp(x_);
    %y_input（1，i）相当于 yk
    s1 = y_input(1,i).*((x_-x_input(1,i+1)).^2).*(delta_h(1,i) + 2.*(x_ - x_input(1,i)))./(delta_h(1,i).^3);
    s2 = y_input(1,i+1).*((x_-x_input(1,i)).^2).*(delta_h(1,i) + 2.*(x_input(1,i+1) - x_))./(delta_h(1,i).^3);
    s3 = m(1,i).*((x_ - x_input(1,i+1)).^2).*(x_ - x_input(1,i))./(delta_h(1,i).^2);
    s4 = m(1,i+1).*((x_ - x_input(1,i)).^2).*(x_ - x_input(1,i+1))./(delta_h(1,i).^2);
    s = s1 + s2 + s3 +s4;
%     disp('s:')
%     disp(s)
    % 判断输入的x属于哪个插值区间，满足则计算对应的f(x)的值
    if x_value>x_input(1,i)&&x_value<x_input(1,i+1)
        s1_ = y_input(1,i)*((x_value-x_input(1,i+1)).^2)*(delta_h(1,i) + 2.*(x_value - x_input(1,i)))./(delta_h(1,i).^3);
        s2_ = y_input(1,i+1)*((x_value-x_input(1,i)).^2)*(delta_h(1,i) + 2.*(x_input(1,i+1) - x_value))./(delta_h(1,i).^3);
        s3_ = m(1,i)*((x_value - x_input(1,i+1)).^2)*(x_value - x_input(1,i))./(delta_h(1,i).^2);
        s4_ = m(1,i+1)*((x_value - x_input(1,i)).^2)*(x_value - x_input(1,i+1))./(delta_h(1,i).^2);
        s_ = s1_ + s2_ + s3_ +s4_;
        % 绘制计算点的结果
        plot(x_value,s_,'ro');
    end
    % 绘制每一段插值函数图像
    plot(x_,s,'b');
    hold on;
end
% 显示网格
plot(x_input,y_input,'*');
disp('xvalue——result')
disp(s_)
axis([0,10,0,6]);
grid on

//main:
x=[1 2 3 4];
y=[2 5 7 9];
m=[-1 2 -6 9];
her_three(x,y,m,3.5)

lagrange插值

%拉格朗日插值命令
function yy=lagrange(x,y,xx)
%Lagrange 插值，求数据（x,y）所表达式的函数在插值点xx处的插值
m=length(x);
n=length(y);
if m~=n
    error('向量x,y长度必须一致');
end
s=0;
for i=1:n
    t=ones(1,length(xx));
    for j=1:n
        if j~=i
            t=t.*(xx-x(j))/(x(i)-x(j));
        end
    end
    s=s+t*y(i);
end

yy=s;
disp(yy)
end

//main
X = [25 40 50 60];
Y = [95 75 63 54];
x=70
yy=lagrange(X,Y,x) %拉格朗日插值

牛顿插值

function p = newton_cha(X,Y,x)
%NEWTON_CHA 此处显示有关此函数的摘要
%   此处显示详细说明

%X,Y是两个行向量，x是要求的插值点的x值

n=length(X);
f=zeros(n,n);

%对cjasj
for k=1:n      
    f(k)=Y(k);
end

%计算各阶插商
for i=2:n  % 差商表从0阶开始；但是矩阵是从1维开始存储！！！！！！
    for k=i:n
        f(k,i)=(f(k,i-1)-f(k-1,i-1))/(X(k)-X(k+1-i));  
    end
end


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

p=f(1,1)+p;
disp(p);

% a = 1980:1:2020;
% p = subs(p,x,a);
% plot(a,p)


% p=0;          
% for k=2:n
%     t=1;
%     for j=1:k-1
%         t=t*(x-X(j));
%     end
%     p=f(k,k)*t+p;
%     disp(p)
% end
% 
% p=f(1,1)+p
fprintf('result is %f\n',p);
end

//main

X=1994:1:2003
Y = [67.052 68.008 69.803 72.024 73.400 72.063 74.669 74.487 74.065 76.777];
x=2010
newton_cha(X,Y,x)

三次样条

X = [0 1 2 3];
Y = [3 5 4 1];
x = 0:0.01:3;
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))*1;
    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)))*1;
    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=[1.000;M;1.000];
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))
         k1=M(i,1)*(X(i+1)-x(t))^3/(6*h(i));
         k2=M(i+1,1)*(x(t)-X(i))^3/(6*h(i));
         k3=(A(i,1)-M(i,1)/6*(h(i))^2)*(X(i+1)-x(t))/h(i);
         k4=(A(i+1,1)-M(i+1,1)/6*(h(i))^2)*(x(t)-X(i))/h(i);
         s(t)=k1+k2+k3+k4; 
         break;
      else
          s(t)=0; 
      end
   end
end
plot(x,s)

牛顿前插公式

function p= Newton_fun(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

//main
 xi=[10 11 12 13 14];
 yi=[2.3026 2.3979 2.4849 2.5649 2.6391];
p= Newton_fun(11.5,xi,yi)