三次样条插值之MATLAB实现

三次样条插值之MATLAB实现

Code

clear all
x = input('输入横坐标，格式:[x1 x2 …… xn]\n')
y = input('输入纵坐标，格式:[y1 y2 …… yn]\n')
n = length(x);

flag = input('请选择边界条件：①已知两端一阶导数值，②已知两端二阶导数值。选择 1 or 2 :');
if flag == 1
    y1_deri = input('x1 的 一阶导数值：');
    yn_deri = input('xn 的 一阶导数值：');
else
    y1_deri = input('x1 的 二阶导数值：');
    yn_deri = input('xn 的 二阶导数值：');
end

for i = 1 : n-1
    h(i) = x(i+1) -  x(i);
end
fprintf('计算 h 结果为:\n');
h

for i = 2 : n-1
    u(i-1) = h(i-1) / (h(i-1) + h(i));
    lamda(i) = h(i) / (h(i-1) + h(i));
end

if flag == 1
    u(n-1) = 1;
    lamda(1) = 1;
else
    u(n-1) = 0;
    lamda(1) = 0;
end

fprintf('计算 μ 结果为: \n');
u
fprintf('计算 λ 结果为：\n');
lamda

for i = 2 : n-1
    d(i) = 6 * ((y(i+1)-y(i))/(x(i+1)-x(i)) - (y(i) - y(i-1))/( x(i) - x(i-1)))/(h(i-1)+h(i));
end

if flag == 1
    d(1) = 6 / h(1) * (( y(2)-y(1) )/( x(2) - x(1) ) - y1_deri);
    d(n) = 6 / h(n-1) * (yn_deri - ((y(n) - y(n-1)) / ( x(n) - x(n-1))));
else
    d(1) = 2 * y1_deri;
    d(n) = 2 * yn_deri;
end

fprintf('计算 d 的结果：\n');
d

matrix1 = zeros(n,n);
for i = 1 : n-1
    matrix1(i,i) = 2;
    matrix1(i,i+1) = lamda(i);
    matrix1(i+1,i) = u(i);
end
matrix1(n,n) = 2;
matrix1;

fprintf('求得 M 结果:\n');
M = matrix1\d'

for i = 1 : n-1
    clear S
    syms t
    k = x(i):0.001:x(i+1);
    fprintf('区间为[ %.3f : %.3f]\n',x(i),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)
    s = M(i)*(x(i+1)-k).^3/ (6*h(i)) + M(i+1)*(k - x(i)).^3/(6*h(i))+(y(i) - M(i)*h(i)^2/6)*(x(i+1) - k)/h(i) + (y(i+1) - M(i+1)*h(i)^2/6)*(k - x(i))/h(i);
    hold on;
    plot(k,s);
end

Result

测试例题：《计算方法》（李庆扬版）P44 例7

>> sanciyangtiao
输入横坐标，格式:[x1 x2 …… xn]
[27.7 28 29 30]

x =

   27.7000   28.0000   29.0000   30.0000

输入纵坐标，格式:[y1 y2 …… yn]
[4.1 4.3 4.1 3.0]

y =

    4.1000    4.3000    4.1000    3.0000

请选择边界条件：①已知两端一阶导数值，②已知两端二阶导数值。选择 1 or 2 :1
x1 的 一阶导数值：3.0
xn 的 一阶导数值：-4.0
计算 h 结果为:

h =

    0.3000    1.0000    1.0000

计算 μ 结果为:

u =

    0.2308    0.5000    1.0000

计算 λ 结果为：

lamda =

    1.0000    0.7692    0.5000

计算 d 的结果：

d =

  -46.6667   -4.0000   -2.7000  -17.4000

求得 M 结果:

M =

  -23.5314
    0.3960
    0.8297
   -9.1149

区间为[ 27.700 : 28.000]

S =

(35650*(t - 28)^3)/2727 - (107*t)/202 + (200*(t - 277/10)^3)/909 + 19317/1010

区间为[ 28.000 : 29.000]

S =

(419*(t - 28)^3)/3030 - (55*t)/202 - (20*(t - 29)^3)/303 + 35929/3030

区间为[ 29.000 : 30.000]

S =

(563*t)/1010 - (4603*(t - 29)^3)/3030 - (419*(t - 30)^3)/3030 - 36977/3030