[数] 三次样条插值的参考程序

本实验项目的要求是让大家通过求解以节点处的一阶导数值为未知量的线性方程组，实现对自然边值条件下的三次样条函数的确定。这个问题看上很麻烦，但是只要把步骤弄清楚，参考我在课堂上给大家的程序，小心仔细的设计那些看似复杂的参数，求解三对角方程即可得到结果。大部分同学完成得很好。我们建议大家在学习完消元法求解线性方程组的知识后，再回头来看一下本程序中线性方程组的求解过程，尝试自己写一个求解三对角方程的程序或函数。
示范程序如下：

代码:

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% spline_student.m
% 本程序用于生成自然边界条件（第2种边界条件）下的样条曲线，并估计函数在某些点处的近似函数值
% 其中x,y为同样长度的向量，代表离散数据。我们要在给定x,y这些数据的前提下，对其进行三次样条插值
% 程序大体分为以下几个步骤：
% 1. 声明(定义)所需的系数矩阵A，右端向量g以及需要的常向量mu,lambda,h
% 2. 为系数矩阵和右端向量赋值
% 3. 求解线性方程组
% 4. 利用线性方程组求解得到的结果，计算样条表达式中的系数C.
% 5. 绘制样条函数图象
% 6. 计算给定点x_nu处的样条函数值，即被逼近的函数在该点处的近似值
%
% note: 本程序是通过求解节点处一阶导数值相关的线性方程组(三对角)实现的
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

% 清空内存，命令行，关闭所有的figure窗口
% clear; clc; close all;
clear; close all;

x = [1 2 4 5 6 7];
y = [1 3 4 2 3 4];

n = length(x);
A = zeros(n,n);
g = zeros(n,1);
mu = zeros(n-2,1);
lambda = zeros(n-2,1);
h = zeros(n-1,1);
h(1) = x(2) - x(1);
for j = 2:n-1
    h(j) = x(j+1) - x(j);   
    lambda(j-1) = h(j)/(h(j) + h(j-1));
    mu(j-1) = 1 - lambda(j-1);
    g(j) = 3*(mu(j-1)*(y(j+1) - y(j))/h(j)+lambda(j-1)*(y(j) - y(j-1))/h(j-1));
end

for j = 2:n-1
    A(j,j-1) = lambda(j-1);
    A(j,j) = 2;
    A(j,j+1) = mu(j-1);
end
A(1,1) = 2;A(1,2)=1;A(n,n)=2;A(n,n-1)=1;
g(1) = 3*(y(2)-y(1))/h(1);
g(n) = 3*(y(n)-y(n-1))/h(n-1);

m = A\g;

C = zeros(n-1,3);
for j = 1:n-1
    C(j,1) = m(j);
    C(j,2) = (3*(y(j+1)-y(j))/h(j)-2*m(j)-m(j+1))/h(j);
    C(j,3) = (m(j+1) + m(j)-2*(y(j+1)-y(j))/h(j))/h(j)^2;
end

h_test = 0.01;
x_test = x(1):h_test:x(n);
m = length(x_test);
y_test = zeros(m,1);
j = 1;
for i = 1:m
    if x_test(i) <= x(j+1)
        y_test(i) = y(j) + C(j,1)*(x_test(i) - x(j)) + C(j,2)*(x_test(i) - x(j))^2 + C(j,3)*(x_test(i) - x(j))^3;
    else
        j = j + 1;
        y_test(i) = y(j) + C(j,1)*(x_test(i) - x(j)) + C(j,2)*(x_test(i) - x(j))^2 + C(j,3)*(x_test(i) - x(j))^3;
    end
end

plot(x,y,'*',x_test,y_test,'r')
legend('插值节点','三次样条',2)
axis auto

x_nu = 3;
for j = 1:n-1
    if x_nu > x(j) && x_nu <= x(j+1)
        break;
    else
        continue;
    end
end
y_nu = y(j) + C(j,1)*(x_nu - x(j)) + C(j,2)*(x_nu - x(j))^2 + C(j,3)*(x_nu - x(j))^3;
y_nu