三次樣條插值matlab實現
%三次樣條差值-matlab通用程序 - zhangxiaolu2015的專欄 - CSDN博客 https://blog.csdn.net/zhangxiaolu2015/article/details/42744823
%【圖文】三次樣條插值算法詳解_百度文庫 https://wenku.baidu.com/view/14423f2e1711cc7931b716ae.html與課堂使用PPT一致。
clc
clear
x=input('請按照格式[x1,x2,x3...]格式輸入y=f(x)函數已知點的橫坐標xi='); %三次樣條差值函數
y=input('請按照格式[y1,y2,y3...]格式輸入y=f(x)函數已知點對應的縱坐標yi=');
x
y
n=size(x,2); %特別注意，matlab中的矩陣編號是從1開始的，而教材上的矩陣編號是從0開始的，即本程序n比PPT上n值大1
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
lambda
以上無論是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
fprintf('邊界條件類型選擇:\n1.已知f(a)和f(b)的二階導數\n2.已知f(a)和f(b)的一階導數\n3.y=f(x)是以T=b-a為周期的周期函數\n');
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
m
1 %三轉角公式 2 3 for k=1:(n-1) 4 5 clear S1 6 7 syms X 8 9 S1=(1-2*(X-x(k))/(-h(k)))*((X-x(k+1))/(h(k)))^2*y(k)+... 10 11 (X-x(k))*((X-x(k+1))/(h(k)))^2*m(k)+... 12 13 (1-2*(X-x(k+1))/(h(k)))*((X-x(k))/(h(k)))^2*y(k+1)+... 14 15 (X-x(k+1))*((X-x(k))/(h(k)))^2*m(k+1); 16 17 fprintf('當%d=<X=<%d時\n',x(k),x(k+1)); 18 19 S=expand(S1) 20 21 end
$$
\begin{array}{l}
{\rm{S(x)}} = {m_k}(X - {x_k}){\left( {\frac{{X - {x_{k + 1}}}}{{{h_k}}}} \right)^2} + \\
{m_{k + 1}}(X - {x_{k + 1}}){\left( {\frac{{X - {x_k}}}{{{h_k}}}} \right)^2} + \\
{y_k}\left( {1 - \frac{{2(X - {x_k})}}{{{-h_k}}}} \right){\left( {\frac{{X - {x_{k + 1}}}}{{{h_k}}}} \right)^2} + \\
{y_{k + 1}}{\left( {\frac{{X - {x_k}}}{{{h_k}}}} \right)^2}\left( {1 - \frac{{2(X - {x_{k + 1}})}}{{{h_k}}}} \right)
\end{array}
$$
