数值分析习题1/牛顿四次插值&自然边界条件下的三次样条插值

题目：牛顿四次插值&自然边界条件下的三次样条插值

牛顿四次插值程序

%%牛顿插值程序.
clc;
format long;                %显示15位
%输入初始数据
X=[0.2 0.4 0.6 0.8 1.0];
Y=[0.98 0.92 0.81 0.64 0.38];
k=1;
N=rand;
for x=0.2:0.08:1.0          %循环使用不同插值点产生函数值数组N(k)
n=max(size(X));             %n  为节点个数
s=1;                        %s  用以迭代每次循环产生（x-xi）的累乘
y=Y(1);            %y 公式第一项
                   %后用以累加新项产生公式
f=Y;               %f  记录差商表的第一列
                   %后用以迭代每次循环产生的差商表新列
                   
for i=1:n-1        %循环产生公式
    f1=f;          %保留新产生的差商表列，以产生下一列
    for j=1:n-i    %循环产生差商表新列
        f(j)=(f1(j+1)-f1(j))/(X(i+j)-X(j));
    end
    fi=f(1);       %fi  公式要用到的i阶差商
    s=s*(x-X(i));
    y=y+s*fi;                     %计算
end

N(k)=y;
k=k+1;
end
x=0.2:0.08:1.0;
plot(x,N,'-r',X,Y,'ob');
xlabel('x');
ylabel('四次NEWTON插值P4（x）');

自然边界条件下的三次样条插值

x=[0.2; 0.4; 0.6; 0.8; 1.0];
y=[0.98; 0.92; 0.81; 0.64; 0.38];
[M,h]=get_M(x,y);
% 求解 M h
% M：三次样条插值函数中的弯矩 
% h：自变量相邻两点差值向量 

% 将 M、h 带入到三次样条插值函数表达式中，从而得出各小区间上的三次样条插值函数 
t = 0.2:0.05:1; 
[~,m] = size(t); 
S=zeros(1,m); 
F=zeros(1,m); 
[n,~] = size(x); 
n=n-1; 
for i=1:n
    tt=(t >=x(i)) & (t <= x(i+1));
    S=( M(i)*((x(i+1) - t).^3./(6*h(i))) ...
        + M(i+1)*((t - x(i)).^3./(6*h(i))) ... 
        + (y(i) - M(i)*(h(i)^2)/6).*((x(i+1) - t)/h(i)) ...
        + (y(i+1) - M(i+1)*(h(i)^2)/6).*((t - x(i))./h(i)) ).* tt;
    for j = 1:m % 此处影响效率，暂如此处理
        if tt(j)
            F(j) = S(j);
        end
    end
    
end 
figure; 
plot(x,y,'b+:'); 
hold on 
plot(t,F,'r');
% hold off 
xlabel('x'),ylabel('y'); 
title('Spline插值'),legend('original points','spline points'); 
grid on; 

get_M函数

% Function：求解自然边界条件下三弯矩方程组对应的矩阵形式的 M A d h
function [M,h] =get_M(x,y) 
[n,~] = size(x); 
n = n-1;

% 自变量相邻两点差值向量
h = diff(x); 
 
% 初始化三对角矩阵A
A = 2 * eye(n-1,n-1);
 
% 求解三对角矩阵A的下对角参数μ
u = zeros(n,1); 
u(n) = 1; 
for i = 1:(n-1) 
    u(i) = h(i) / (h(i) + h(i+1)); 
end 
 
% 求解三对角矩阵A的上对角参数λ
lambda = zeros(n,1); 
lambda(1) = 1; 
for i = 2:n 
    lambda(i) = h(i) / (h(i-1) + h(i)); 
end 


% 求解矩阵方程右端向量d
d = zeros(n-1,1); 
for i = 2:n 
    d(i-1) = 6 / (x(i+1) - x(i-1)) * ((y(i+1) - y(i))/(x(i+1) - x(i)) - (y(i) - y(i-1))/(x(i) - x(i-1))); 
end 


u1=u(2:n-1);
lambda1=(2:n-1);
% 求解三对角矩阵A
A = A + diag(u1,-1) + diag(lambda1,1); 
% 矩阵方程左乘A的逆矩阵求解出M 
M1= A^(-1)*d;
M(1)=0;
M(n+1)=0;
for k=2:n
    M(k)=M1(k-1);
end

end