1.三次样条插值函数
%%三次样条插值
%%bc为boundary conditions(边界条件),当已知两端点的一阶导数值时为-1,当已知两端的二阶导数时为0,当函数为周期函数时为1
%%X为节点值,Y为函数表达式(attribute=0)或者具体值(attribute=1)
function CSI = Cubic_spline_interpolation(X,Y,precision,attribute,bc)
[m,n] = size(X);a = min(X);b = max(X);n = n-1;
X = sort(X);
%%画已知函数或已知点图像
if attribute == 0
F = subs(Y,X);
t = a:(b-a)/precision:b;
Y_real = subs(Y,t);
pic = figure;
set(pic,'color','w');
plot(X,F,'r*',t,Y_real,'b');
grid on
xlabel('x shaft');ylabel('y shaft');
title('三次样条插值');
hold on
elseif attribute ==1
F = Y;
pic = figure;
set(pic,'color','w');
plot(X,F,'r*');
grid on
xlabel('x shaft');ylabel('y shaft');
title('三次样条插值');
hold on
end
if bc == -1
left_endpoint = input('输入左端点的一阶导数值:');
right_endpoint = input('输入右端点的一阶导数值:');
for i = 1:n
X_one_interval{i} = [X(i),X(i+1)];%%节点构造区间
F_one_interval{i} = [F(i),F(i+1)];%%构造节点值区间
F_all_bad(i) = The_first_order_All_bad(F_one_interval{i},X_one_interval{i},attribute);%%节点区间的一阶均差
h(i) = X(i+1)-X(i);
end
%%样条插值算法
for i = 1:n-1
mu(i) = h(i)/(h(i)+h(i+1));
end
mu(n) = 1;
lambda(1) = 1;
for i = 2:n
lambda(i) = h(i)/(h(i-1)+h(i));
end
d(1)=6*(F_all_bad(1)-left_endpoint)/h(1);
for i = 2:n
d(i) = 6*(F_all_bad(i)-F_all_bad(i-1))/(h(i-1)+h(i));
end
d(n+1) = 6*(right_endpoint-F_all_bad(n))/h(n);
A = zeros(n+1,n+1);
for i = 1:n+1
for j = 1:n+1
if i == j
A(i,j) = 2;
elseif (i-j) == 1
A(i,j) = mu(j);
elseif (j-i) == 1
A(i,j) = lambda(i);
end
end
end
%%解样条算法中的线性方程组得出解
disp('样条函数初始系数:');
M = vpa((A\d')',6)
%%输出节点区间及对应的样条函数
syms x;
for i = 1:n
CSI{1,i} = X_one_interval{i};
CSI{2,i} = vpa((X(i+1)-x)^3*M(i)/(6*h(i))+(x-X(i))^3*M(i+1)/(6*h(i))+(F(i)-M(i)*h(i)^2/6)*(X(i+1)-x)/h(i)+(F(i+1)-M(i+1)*h(i)^2/6)*(x-X(i))/h(i),5);
end
%%画样条函数的图像
for i = 1:n
t_pic{i} = X(i):(X(i+1)-X(i))/precision:X(i+1);
T_pic{i} = subs(CSI{2,i},t_pic{i});
end
for i = 1:n
t_Pic(1+(precision+1)*(i-1):(precision+1)*i) = t_pic{i};
T_Pic(1+(precision+1)*(i-1):(precision+1)*i) = T_pic{i};
end
plot(t_Pic,T_Pic,'g');
if attribute ==0
legend('F:数据点','Y_real:真实图像','T_Pic:样条插值拟合函数');
elseif attribute == 1
legend('F:数据点','T_Pic:样条插值拟合函数');
end
elseif bc == 0
left_endpoint = input('输入左端点的二阶导数值:');
right_endpoint = input('输入右端点的二阶导数值:');
for i = 1:n
X_one_interval{i} = [X(i),X(i+1)];
F_one_interval{i} = [F(i),F(i+1)];
F_all_bad(i) = The_first_order_All_bad(F_one_interval{i},X_one_interval{i},attribute);
h(i) = X(i+1)-X(i);
end
for i = 1:n-1
mu(i) = h(i)/(h(i)+h(i+1));
end
mu(n) = 0;
lambda(1) = 0;
for i = 2:n
lambda(i) = h(i)/(h(i-1)+h(i));
end
d(1)=2*left_endpoint;
for i = 2:n
d(i) = 6*(F_all_bad(i)-F_all_bad(i-1))/(h(i-1)+h(i));
end
d(n+1) = 2*right_endpoint;
A = zeros(n+1,n+1);
for i = 1:n+1
for j = 1:n+1
if i == j
A(i,j) = 2;
elseif (i-j) == 1
A(i,j) = mu(j);
elseif (j-i) == 1
A(i,j) = lambda(i);
end
end
end
disp('样条函数初始系数:');
M = vpa((A\d')',6)
syms x;
for i = 1:n
CSI{1,i} = X_one_interval{i};
CSI{2,i} = vpa((X(i+1)-x)^3*M(i)/(6*h(i))+(x-X(i))^3*M(i+1)/(6*h(i))+(F(i)-M(i)*h(i)^2/6)*(X(i+1)-x)/h(i)+(F(i+1)-M(i+1)*h(i)^2/6)*(x-X(i))/h(i),5);
end
for i = 1:n
t_pic{i} = X(i):(X(i+1)-X(i))/precision:X(i+1);
T_pic{i} = subs(CSI{2,i},t_pic{i});
end
for i = 1:n
t_Pic(1+(precision+1)*(i-1):(precision+1)*i) = t_pic{i};
T_Pic(1+(precision+1)*(i-1):(precision+1)*i) = T_pic{i};
end
plot(t_Pic,T_Pic,'g');
if attribute ==0
legend('F:数据点','Y_real:真实图像','T_Pic:样条插值拟合函数');
elseif attribute == 1
legend('F:数据点','T_Pic:样条插值拟合函数');
end
end
end
2.一阶均差函数
%%一阶均差
function FAb = The_first_order_All_bad(f,X,attribute)
[m,n] = size(X);
if attribute == 0
F = subs(f,X);
FAb = (F(n)-F(1))/(X(n)-X(1));
elseif attribute == 1
FAb = (f(n)-f(1))/(X(n)-X(1));
end
end
3.样条插值拟合值函数
%%三次样条插值拟合值
%%用法同三次样条函数
function CSIV = Cubic_spline_interpolation_value(X,Y,precision,attribute,bc,x_value)
[m,n] = size(X);a = min(X);b = max(X);n = n-1;
X = sort(X);
if attribute == 0
F = subs(Y,X);
elseif attribute ==1
F = Y;
end
if bc == -1
left_endpoint = input('输入左端点的一阶导数值:');
right_endpoint = input('输入右端点的一阶导数值:');
for i = 1:n
X_one_interval{i} = [X(i),X(i+1)];
F_one_interval{i} = [F(i),F(i+1)];
F_all_bad(i) = The_first_order_All_bad(F_one_interval{i},X_one_interval{i},attribute);
h(i) = X(i+1)-X(i);
end
for i = 1:n-1
mu(i) = h(i)/(h(i)+h(i+1));
end
mu(n) = 1;
lambda(1) = 1;
for i = 2:n
lambda(i) = h(i)/(h(i-1)+h(i));
end
d(1)=6*(F_all_bad(1)-left_endpoint)/h(1);
for i = 2:n
d(i) = 6*(F_all_bad(i)-F_all_bad(i-1))/(h(i-1)+h(i));
end
d(n+1) = 6*(right_endpoint-F_all_bad(n))/h(n);
A = zeros(n+1,n+1);
for i = 1:n+1
for j = 1:n+1
if i == j
A(i,j) = 2;
elseif (i-j) == 1
A(i,j) = mu(j);
elseif (j-i) == 1
A(i,j) = lambda(i);
end
end
end
M = (A\d')';
syms x;
for i = 1:n
S{1,i} = X_one_interval{i};
S{2,i} = vpa((X(i+1)-x)^3*M(i)/(6*h(i))+(x-X(i))^3*M(i+1)/(6*h(i))+(F(i)-M(i)*h(i)^2/6)*(X(i+1)-x)/h(i)+(F(i+1)-M(i+1)*h(i)^2/6)*(x-X(i))/h(i),5);
end
for i =1:n
if x_value >= X(i) && x_value <= X(i+1);
s = i;
end
end
CSIV{1,1} = '拟合值';
CSIV{2,1} = vpa(subs(S{2,s},x_value),4);
CSIV{1,2} = '实际值';
CSIV{2,2} = vpa(subs(Y,x_value),4);
CSIV{1,3} = '误差';
CSIV{2,3} = abs(CSIV{2,1}-CSIV{2,2});
elseif bc == 0
left_endpoint = input('输入左端点的二阶导数值:');
right_endpoint = input('输入右端点的二阶导数值:');
for i = 1:n
X_one_interval{i} = [X(i),X(i+1)];
F_one_interval{i} = [F(i),F(i+1)];
F_all_bad(i) = The_first_order_All_bad(F_one_interval{i},X_one_interval{i},attribute);
h(i) = X(i+1)-X(i);
end
for i = 1:n-1
mu(i) = h(i)/(h(i)+h(i+1));
end
mu(n) = 0;
lambda(1) = 0;
for i = 2:n
lambda(i) = h(i)/(h(i-1)+h(i));
end
d(1)=2*left_endpoint;
for i = 2:n
d(i) = 6*(F_all_bad(i)-F_all_bad(i-1))/(h(i-1)+h(i));
end
d(n+1) = 2*right_endpoint;
A = zeros(n+1,n+1);
for i = 1:n+1
for j = 1:n+1
if i == j
A(i,j) = 2;
elseif (i-j) == 1
A(i,j) = mu(j);
elseif (j-i) == 1
A(i,j) = lambda(i);
end
end
end
M = (A\d')';
syms x;
for i = 1:n
S{1,i} = X_one_interval{i};
S{2,i} = vpa((X(i+1)-x)^3*M(i)/(6*h(i))+(x-X(i))^3*M(i+1)/(6*h(i))+(F(i)-M(i)*h(i)^2/6)*(X(i+1)-x)/h(i)+(F(i+1)-M(i+1)*h(i)^2/6)*(x-X(i))/h(i),5);
end
CSIV{1,1} = '拟合值';
CSIV{2,1} = vpa(subs(S{2,s},x_value),4);
end
4.例子
clear all
clc
syms x;
Y=sin(x)/(1+x^2)+exp(-x^2);
X=-5:1:5;
precision=500;
attribute=0;
bc=-1;
%%三次样条插值
Cubic_spline_interpolation(X,Y,precision,attribute,bc)
结果为
输入左端点的一阶导数值:subs(diff(Y),-5)
输入右端点的一阶导数值:subs(diff(Y),5)
样条函数初始系数:
M =
[ -0.0869325, 0.0232929, -0.00623918, 0.00166378, -0.000415945, 1.32899e-12, 0.000415945, -0.00166378, 0.00623918, -0.0232929, 0.0869325]
S =
2×10 cell 数组
列 1 至 4
{1×2 double} {1×2 double} {1×2 double} {1×2 double}
{1×1 sym } {1×1 sym } {1×1 sym } {1×1 sym }
列 5 至 8
{1×2 double} {1×2 double} {1×2 double} {1×2 double}
{1×1 sym } {1×1 sym } {1×1 sym } {1×1 sym }
列 9 至 10
{1×2 double} {1×2 double}
{1×1 sym } {1×1 sym }
>> S{2,:}
ans =
0.014489*(x + 4.0)^3 - 0.010735*x + 0.0038822*(x + 5.0)^3 - 0.0023031
ans =
- 0.053584*x - 0.0010399*(x + 4.0)^3 - 0.0038822*(x + 3.0)^3 - 0.1737
ans =
0.0010399*(x + 2.0)^3 - 0.15087*x + 0.0002773*(x + 3.0)^3 - 0.46557
ans =
0.11103*x - 0.0002773*(x + 1.0)^3 - 0.000069324*(x + 2.0)^3 + 0.058248
ans =
1.0528*x + 2.215e-13*(x + 1.0)^3 + 0.000069324*x^3 + 1.0
ans =
0.000069324*x^3 - 2.215e-13*(x - 1.0)^3 - 0.21145*x + 1.0
ans =
1.3766 - 0.0002773*(x - 1.0)^3 - 0.000069324*(x - 2.0)^3 - 0.58809*x
ans =
0.0010399*(x - 2.0)^3 - 0.18726*x + 0.0002773*(x - 3.0)^3 + 0.57497
ans =
0.17469 - 0.0010399*(x - 4.0)^3 - 0.0038822*(x - 3.0)^3 - 0.053831*x
ans =
0.014489*(x - 4.0)^3 - 0.010735*x + 0.0038822*(x - 5.0)^3 + 0.0023042
>>
函数图像为