三次样条插值matlab实现

%三次样条差值-matlab通用程序 - zhangxiaolu2015的专栏 - CSDN博客 https://blog.csdn.net/zha
%【图文】三次样条插值算法详解_百度文库 https://wenku.baidu.com/view/14423f2e1711cc7931b716
clc
clear
x=input('请按照格式[x1,x2,x3...]格式输入y=f(x)函数已知点的横坐标xi='); %三次样条差值函数
y=input('请按照格式[y1,y2,y3...]格式输入y=f(x)函数已知点对应的纵坐标yi=');
x

x = 1x4 double
1 2 4 5

y

y = 1x4 double
1 3 4 2

n=size(x,2); %特别注意，matlab中的矩阵编号是从1开始的，而教材上的矩阵编号是从0开始的，即本程序
for k=2:n %计算h(i)
h(k-1)=x(k)-x(k-1);
end
for k=1:(n-2) %计算μ和λ
mu(k)=h(k)/(h(k)+h(k+1));
lambda(k)=1-mu(k);
end
mu

mu = 1x2 double
0.3333 0.6667

lambda

lambda = 1x2 double
0.6667 0.3333

以上无论是M还是m关系式矩阵通用。

for k=1:(n-2)
  g(k)=3*(lambda(k)*(y(k+1)-y(k))/h(k)+mu(k)*(y(k+2)-y(k+1))/h(k+1));      %计算g(1)到g(n-2)
end
g

g =

-1.288728000000000 -2.093712750000000 -3.177727125000001

fprintf('边界条件类型选择:\n1.已知f(a)和f(b)的二阶导数\n2.已知f(a)和f(b)的一阶导数\n');

边界条件类型选择:

1.已知f(a)和f(b)的二阶导数

2.已知f(a)和f(b)的一阶导数

3.y=f(x)是以T=b-a为周期的周期函数

in=input('请输入对应序号：');
if in==1
    in
    M(1)=input('请输入f(a)的二阶导数值：');
    M(n)=input('请输入f(b)的二阶导数值：');
    M(1)
    M(n)
    A=zeros(n,n); %构造追赶法所需的A和b
for k=2:(n-1)
    A(k,k)=2;
    A(k,k+1)=mu(k-1);
    A(k,k-1)=lambda(k-1);
end
    A(1,1)=2;
    A(1,2)=1;
    A(end,end)=2;
    A(end,end-1)=1;
    A
    b=zeros(n,1);
for k=2:(n-1)
    b(k,1)=g(k-1);
end
    b(1,1)=3*((y(2)-y(1))/h(1)-2*h(1)*M(1));
    b(n,1)=3*((y(n)-y(n-1))/h(n-1)+2*h(n-1)*M(n));
    b
    b=b';
    m=zhuigan(A,b); %利用追赶法求解成功，这里的参数b形式应为行向量而非列向量
elseif in==2
    y0=input('请输入f(a)的一阶导数值：');
    yn=input('请输入f(b)的一阶导数值：');
    A=zeros(n-2,n-2); %构造追赶法所需的A和b
for k=2:(n-3)
    A(k,k)=2;
    A(k,k+1)=mu(k);
    A(k,k-1)=lambda(k);
end
    A(1,1)=2;
    A(1,2)=mu(1);
    A(end,end)=2;
    A(end,end-1)=lambda(n-2);
    b=zeros(n-2,1);
for k=2:(n-3)
    b(k,1)=g(k);
end
    b(1,1)=g(1)-lambda(1)*y0;
    b(end,1)=g(n-2)-mu(n-2)*yn;
    b=b';
    m=zhuigan(A,b);%利用追赶法求解
    m(1)
    m(2)
%这里解出m(1)至m(n-2),为能代入带一阶导数的分段三次埃米尔特插值多项式，要对m进行调整
for k=(n-2):-1:1
    m(k+1)=m(k);
end
    m(1)=y0;
    m(n)=yn;
elseif in==3
    A=zeros(n,n); %构造追赶法所需的A和b
for k=2:(n-1)
    A(k,k)=2;
    A(k,k+1)=mu(k-1);
    A(k,k-1)=lambda(k-1);
end
    A(1,1)=2;
    A(1,2)=mu(1);
    A(1,end)=lambda(1);
    A(end,end)=2;
    A(end,end-1)=lambda(n-1);
    A(end,1)=mu(n-1);
    b=zeros(n-1,1);
for k=1:(n-1)
    b(k,1)=d(k+1);
end
    N=LU_fenjieqiuxianxingfangcheng(A,b); %利用LU分解求解线性方程组
for k=1:(n-1)
    M(k+1)=N(k,1);
end
    M(1)=M(n);
else
    fprintf('您输入的序号不正确');
end

ans = 0

A = 4x4 double

​ 2.0000 1.0000 0 0
​ 0.6667 2.0000 0.3333 0
​ 0 0.3333 2.0000 0.6667
​ 0 0 1.0000 2.0000
b = 4x1 double

​ 6.0000
​ 4.5000
​ -3.5000
​ -6.0000
c = 1x3 double

​ 0.6667 0.3333 1.0000
a = 1x4 double

2 2 2 2
b = 1x3 double

1.0000 0.3333 0.6667

m

m = 1x4 double
2.1250 1.7500 -1.2500 -2.3750

%三转角公式
for k=1:(n-1)
clear S1
syms X
S1=(1-2*(X-x(k))/(-h(k)))*((X-x(k+1))/(h(k)))^2*y(k)+...
(X-x(k))*((X-x(k+1))/(h(k)))^2*m(k)+...
(1-2*(X-x(k+1))/(h(k)))*((X-x(k))/(h(k)))^2*y(k+1)+...
(X-x(k+1))*((X-x(k))/(h(k)))^2*m(k+1);
fprintf('当%d=<X=<%d时\n',x(k),x(k+1));
S=expand(S1)
end

当1=<X=<2时

S =
当2=<X=<4时

S =
当4=<X=<5时

S =