MATLAB 样条曲线拟合代码解析及优化

clear, clc  % 清空命令行窗口和工作区变量
X=-1:0.2:1;  % 定义X向量
Y=1./(25*X.^2+1);  % 定义Y向量
dx0= 0.0739644970414201;  % 定义左端点一阶导数
dxn=-0.0739644970414201;  % 定义右端点一阶导数
S=csfit(X,Y,dx0,dxn);  % 调用csfit函数拟合样条曲线
x1=-1:0.01:-0.5;y1=polyval(S(1,:),x1-X(1));  % 计算第一段样条曲线的x和y值
x2=-0.5:0.01:0;y2=polyval(S(2,:),x2-X(2));  % 计算第二段样条曲线的x和y值
x3=0:0.01:0.5; y3=polyval(S(3,:),x3-X(3));  % 计算第三段样条曲线的x和y值
x4=0.5:0.01:1;y4=polyval(S(4,:),x4-X(4));  % 计算第四段样条曲线的x和y值
plot(x1,y1,x2,y2,x3,y3,x4, y4, X,Y,'.')  % 绘制样条曲线和原始数据点
function S=csfit (X,Y,dx0,dxn)  % 定义样条拟合函数
N=length (X)-1;  % 计算样本点个数
H=diff(X);  % 计算相邻样本点之间的距离
D=diff (Y)./H;  % 计算相邻样本点的斜率
A=H(2:N-1);  % 定义矩阵A
B=2*(H(1:N-1)+H(2:N));  % 定义矩阵B
C=H(2:N);  % 定义矩阵C
U=6*diff(D);  % 定义矩阵U
B (1)=B(1)-H(1)/2;  % 调整B矩阵第一个元素
U(1)=U(1)-3*(D(1));  % 调整U矩阵第一个元素
B(N-1)=B(N-1)-H(N)/2;  % 调整B矩阵最后一个元素
U(N-1)=U(N-1)-3*(-D (N));  % 调整U矩阵最后一个元素
for k=2:N-1  % 循环计算A、B、U矩阵
    temp=A (k-1)/B(k-1);
    B(k)=B(k)-temp*C(k-1);
    U(k)=U(k)-temp*U(k-1);
end
M(N)=U(N-1)/B(N-1);  % 计算M矩阵最后一个元素
for k=N-2:-1:1  % 循环计算M矩阵
    M(k+1)=(U(k)-C(k)*M(k+2))/B(k);
end
M(1) =3*(D (1)-dx0)/H(1)-M(2)/2;  % 计算M矩阵的第一个元素
M(N+1)=3* (dxn-D(N))/H(N)-M(N)/2;  % 计算M矩阵的最后一个元素
for k=0:N-1  % 循环计算S矩阵
    S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1));
    S(k+1,2)=M(k+1)/2;
    S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6;
    S(k+1,4)=Y(k+1);
end
end  % 结束函数定义

代码逐句解析内容：

clear, clc % 清空命令行窗口和工作区变量
X=-1:0.2:1; % 定义X向量，从-1到1，步长为0.2
Y=1./(25*X.^2+1); % 定义Y向量，根据X向量计算对应的Y值
dx0= 0.0739644970414201; % 定义左端点一阶导数
dxn=-0.0739644970414201; % 定义右端点一阶导数
S=csfit(X,Y,dx0,dxn); % 调用csfit函数拟合样条曲线，并将结果存储到变量S中
x1=-1:0.01:-0.5;y1=polyval(S(1,:),x1-X(1)); % 计算第一段样条曲线的x和y值，使用polyval函数进行多项式插值
x2=-0.5:0.01:0;y2=polyval(S(2,:),x2-X(2)); % 计算第二段样条曲线的x和y值
x3=0:0.01:0.5; y3=polyval(S(3,:),x3-X(3)); % 计算第三段样条曲线的x和y值
x4=0.5:0.01:1;y4=polyval(S(4,:),x4-X(4)); % 计算第四段样条曲线的x和y值
plot(x1,y1,x2,y2,x3,y3,x4, y4, X,Y,'.') % 绘制样条曲线和原始数据点

csfit 函数解析：

function S=csfit (X,Y,dx0,dxn) % 定义样条拟合函数
N=length (X)-1; % 计算样本点个数，减1是因为需要计算N-1个样条段
H=diff(X); % 计算相邻样本点之间的距离
D=diff (Y)./H; % 计算相邻样本点的斜率
A=H(2:N-1); % 定义矩阵A，用于存储相邻样本点距离，除了第一个和最后一个元素
B=2*(H(1:N-1)+H(2:N)); % 定义矩阵B，用于存储相邻样本点距离的累加
C=H(2:N); % 定义矩阵C，用于存储相邻样本点距离，除了第一个元素
U=6*diff(D); % 定义矩阵U，用于存储相邻样本点斜率的差值
B (1)=B(1)-H(1)/2; % 调整B矩阵第一个元素，考虑边界条件
U(1)=U(1)-3*(D(1)); % 调整U矩阵第一个元素，考虑边界条件
B(N-1)=B(N-1)-H(N)/2; % 调整B矩阵最后一个元素，考虑边界条件
U(N-1)=U(N-1)-3*(-D (N)); % 调整U矩阵最后一个元素，考虑边界条件
for k=2:N-1 % 循环计算A、B、U矩阵，使用高斯消元法求解线性方程组 temp=A (k-1)/B(k-1); B(k)=B(k)-temp*C(k-1); U(k)=U(k)-temp*U(k-1); end
M(N)=U(N-1)/B(N-1); % 计算M矩阵最后一个元素，使用回代法求解线性方程组
for k=N-2:-1:1 % 循环计算M矩阵 M(k+1)=(U(k)-C(k)*M(k+2))/B(k); end
M(1) =3*(D (1)-dx0)/H(1)-M(2)/2; % 计算M矩阵的第一个元素，考虑边界条件
M(N+1)=3* (dxn-D(N))/H(N)-M(N)/2; % 计算M矩阵的最后一个元素，考虑边界条件
for k=0:N-1 % 循环计算S矩阵，S矩阵存储样条曲线的系数 S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1)); S(k+1,2)=M(k+1)/2; S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6; S(k+1,4)=Y(k+1); end
end % 结束函数定义

优化建议：

可以将代码中的 dx0 和 dxn 作为函数参数，以便可以根据需要改变边界条件
可以将 csfit 函数的代码封装成一个单独的 M 文件，方便调用和维护
可以使用 MATLAB 的 spline 函数进行样条曲线拟合，简化代码编写

总结：

本文详细解析了 MATLAB 样条曲线拟合代码，并提供了一些优化建议。希望能帮助读者更好地理解和应用样条曲线拟合方法。

相关链接：