- 📚 相关推荐文章
- 空间插值 推荐
- 拉格朗日插值公式 推荐
- 样条插值 推荐
- matlab编写拉格朗日插值代码函数 推荐
- 插值 推荐
样条插值
-
2023年3月5日发(作者：蓝牙硬件).
;.
三次样条插值
1.算法原理
由于在许多实际问题中，要求函数的二阶导数连续，人们便提出了三次样条插值函数，
三次样条插值函数是由分段三次函数拼接而成的，在连接点处二阶导数连续。
设S(x)在节点i
x
处的二阶导数
),1,0()(''niMxS
ii
，
，其中i
M
为待定参数。由
S（x）是分段三次多项式可知，
)(''xS
是分段线性函数，
)(''xS
在子区间
ii
xx,
1上可以表
示为
iii
i
i
i
i
i
i
ii
i
i
ii
ixxxM
h
xx
M
h
xx
M
xx
xx
M
xx
xx
xS
1
1
1
1
1
1
1
,)(''
其中1
iii
xxh
，对S（x）两端积分两次得
iiiiiii
i
i
i
i
ixxxxxcxxbM
h
xx
M
h
xx
xS
11
1
1
1
3
),(
6
)(
6
)(
)(
其中i
b
和i
c
为积分常数。由插值条件
iiii
yxSyxS
)(,
11得
iiii
i
iiii
iyhcM
h
yhbM
h
6
,
6
2
11
2
由此解得
ii
i
iiii
i
ii
hM
h
ychM
h
yb/
6
,/
6
2
1
2
1
代入得：
ii
i
i
i
i
i
i
i
i
i
ii
i
i
i
i
ixxx
h
xx
M
h
y
h
xx
M
h
yM
h
xx
M
h
xx
xS
1
1
2
1
2
1
3
1
1
3
,
6666
求导得：
iii
ii
i
ii
i
i
i
i
i
ixxxh
MM
h
yy
M
h
xx
M
h
xx
xS
1
11
2
1
1
2
,
622
'
令
i
xx
得xS在
i
x处的左导数
i
ii
i
i
i
i
i
ii
i
ii
i
i
h
yy
M
h
M
h
h
MM
h
yy
M
h
xS1
1
11
3662
'
.
;.
又令
1
i
xx得xS在
1i
x
处右导数
i
ii
i
i
i
i
i
ii
i
ii
i
i
ih
yy
M
h
M
h
h
MM
h
yy
M
h
xS1
1
11
116362
'
，
从而有
1
1
1
11
63
)('
i
ii
i
i
i
i
h
yy
M
h
M
h
xS，由xS在节点
i
x处一阶导数的连续性知
1,2,1),(')('
nixSxS
ii
，，即
1,2,1
636
1
1
1
1
11
1
ni
h
yy
h
yy
M
h
M
hh
M
h
i
ii
i
ii
i
i
i
ii
i
i，，
两端同乘
1
6
ii
hh
得)(
6
21
1
1
1
1
1
1
1
1i
ii
i
ii
ii
i
ii
i
ii
ii
i
h
yy
h
yy
hh
M
hh
h
MM
hh
h
，
记
i
ii
i
i
ii
i
ihh
h
hh
h
1,
1
1
1
，
1,,2,1],,,[6
6
11
1
1
1
1
nixxxf
h
yy
h
yy
hh
d
iii
i
ii
i
ii
ii
i
，则关于
i
M
的方
程组写成
1,2,1,2
11
nidMMM
iiiii
，。
三种边界条件的三弯矩方程：
（1）第一种边界条件：已知bfaf'',''。取)('',''
0
bfMafM
n
，这时方程
组减少了两个未知量，变成只含n-1个未知量1,2,1niM
i
，
的n-1个方程的方程
组
2,3,2,
2
2
2
11121
11
02
ni
MdMM
dMMM
MdMM
nnnnnn
iiiii
iii
，
，用矩阵表示为
nnn
n
n
n
n
nn
Md
d
d
Md
M
M
M
M
11
2
2
011
1
2
2
1
1
22
22
1
2
2
2
2
可用追赶法求解出
)1,2,1(niM
i
，
后，即得三条样条插值函数。
.
;.
（2）第二种边界条件，已知)(')('bfaf，，记
)(''),(''
0
bfyafy
n
，则有
')(',')('
00nn
yxSyxS
，得
'
36
'
63
-1
10
1
01
1
1
0
1
n
n
nn
n
n
n
ny
h
yy
M
h
M
h
y
h
yy
M
h
M
h
，，即
nnn
dMMdMM
2,2
1010
，其中],,[6],,,[6
12100nnnn
xxxfdxxxfd
，得到方
程组
1,2,1,
2
2
2
1
11
010
ni
dMM
dMMM
dMM
nnn
iiiiii
，
，用矩阵表示为
n
n
n
n
nn
d
d
d
d
M
M
M
M
1
1
0
1
1
0
11
11
21
2
2
12
，该方程组的系数矩阵是严格三对角占
优矩阵，可用追赶法求解。
（3）第三种边界条件：周期型边界条件。已知xfy是以
0
xxabT
n
为周
期的周期函数，则由周期性可知，
11110110
,,,hhMMMMyyyy
nnnnn
，
，这
时将点
n
x看成是内节点，则有
nnnnnn
dMMM
11
2，也即
nnnnn
dMMM
111
2，其中
n
nnn
n
nh
yy
h
yy
hh
d1
1
1
1
6
，方程组第1个方
程为：
12111
2dMMM
n
，所以方程组为
1,3,2,
2
2
2
11
11
11211
ni
dMMM
dMMM
dMMM
nnnnn
iiiiii
n
，
，用矩阵表示为
.
;.
n
n
n
n
nn
nn
d
d
d
d
M
M
M
M
1
2
1
1
2
1
11
22
11
2
2
2
2
，显然系数矩阵为严格对角占优矩阵，
可用LU分解法求解。
.
;.
2.程序框图
计算步长
得出三弯矩方程组
解矩阵方程
计算三次样条曲线
的系数
用插值函数计算出
各点的插值计算值
3.源程序
functionx=mchase(A,d)
%追赶法
n=length(d);
u=zeros(n,1);u(1)=A(1,1);
fork=2:n
l(k)=A(k,k-1)/u(k-1);
u(k)=A(k,k)-l(k)*A(k-1,k);
end
y=zeros(n,1);y(1)=d(1);
fori=2:n
y(i)=d(i)-l(i)*y(i-1);
end
x=zeros(n,1);
.
;.
x(n)=y(n)/u(n);
fori=n-1:-1:1
x(i)=(y(i)-A(i,i+1)*x(i+1))/u(i);
end
x
end
functionT=mspline1(x0,y0,f21,f22,xx)
%三次样条插值函数第一种边界条件
%x0、y0分别为节点的横坐标和纵坐标；
%f21为左端点的二阶导数值,f22为右端点的二阶导数值；xx为由插值点组成的向量
n=length(x0)-1;%计算小区间数
fori=1:n
h(i)=x0(i+1)-x0(i);
end
fori=1:n-1
mu(i)=h(i)/(h(i)+h(i+1));
lamda(i)=1-mu(i);
d(i)=6*((y0(i+2)-y0(i+1))/h(i+1)-(y0(i+1)-y0(i))/h(i))/(h(i)+h(i+1));
end
A=zeros(n-1);
fori=1:n-2
A(i+1,i)=mu(i+1);%次下对角线
A(i,i+1)=lamda(i);%次上对角线
A(i,i)=2;%主对角线
end
A(n-1,n-1)=2;
dd=zeros(n-1,1);%右端列向量
fori=2:n-2
dd(i)=d(i);
end
dd(1)=d(1)-mu(1)*f21;dd(n-1)=d(n-1)-lamda(n-1)*f22;
M=mchase(A,dd);%追赶法求解M值
h
mu
lamda
A
dd
M=[f21,M',f22]'
t=sym('t');
a=zeros(n,1);b=zeros(n,1);c=zeros(n,1);e=zeros(n,1);
fori=1:n
a(i)=M(i)./(6*h(i));
b(i)=M(i+1)./(6*h(i));
.
;.
W1(i)=b(i)-a(i);
W2(i)=3*(a(i).*x0(i+1)-b(i).*x0(i));
c(i)=y0(i)./h(i)-h(i).*M(i)/6;
e(i)=y0(i+1)./h(i)-h(i).*M(i+1)/6;
W3(i)=3*b(i).*x0(i).^2-3*a(i).*x0(i+1).^2+e(i)-c(i);
W4(i)=a(i).*x0(i+1).^3-b(i).*x0(i).^3+c(i).*x0(i+1)-e(i).*x0(i);
Si(t)=W1(i).*t^3+W2(i).*t^2+W3(i).*t+W4(i)%每个小区间的三次样条差值函数表达式
end
m=length(xx);T=zeros(m,1);
fork=1:m
forj=1:n
if((xx(k)>=x0(j))&(xx(k) T(k)=W1(j).*(xx(k).^3)+W2(j).*(xx(k).^2)+W3(j).*xx(k)+W4(j); end end end T End functionT=mspline2(x0,y0,f11,f12,xx) %三次样条插值函数第二种边界条件 %x0、y0分别为节点的横坐标和纵坐标； %f11为左端点的二阶导数值,f12为右端点的二阶导数值；xx为由插值点组成的向量 n=length(x0)-1;%计算小区间数 fori=1:n h(i)=x0(i+1)-x0(i); end fori=1:n-1 mu(i)=h(i)/(h(i)+h(i+1)); lamda(i)=1-mu(i); d(i)=6*((y0(i+2)-y0(i+1))/h(i+1)-(y0(i+1)-y0(i))/h(i))/(h(i)+h(i+1)); end A=zeros(n+1);%系数矩阵 dd=zeros(n+1,1);%右端列向量 fork=2:n A(k,k)=2;%主对角线元素 A(k,k-1)=mu(k-1);%次下对角线元素 A(k,k+1)=lamda(k-1);%次上对角线元素 end A(1,1)=2;A(1,2)=1;A(n+1,n+1)=2;A(n+1,n)=1; dd(1)=6*((y0(2)-y0(1))/h(1)-f11)/h(1); dd(n+1)=6*(f12-(y0(n+1)-y0(n))/h(n))/h(n); fork=2:n . ;. dd(k)=d(k-1); end M=mchase(A,dd);%追赶法求解M值 h mu lamda A dd M t=sym('t'); a=zeros(n,1);b=zeros(n,1);c=zeros(n,1);e=zeros(n,1); fori=1:n a(i)=M(i)./(6*h(i)); b(i)=M(i+1)./(6*h(i)); W1(i)=b(i)-a(i); W2(i)=3*(a(i).*x0(i+1)-b(i).*x0(i)); c(i)=y0(i)./h(i)-h(i).*M(i)/6; e(i)=y0(i+1)./h(i)-h(i).*M(i+1)/6; W3(i)=3*b(i).*x0(i).^2-3*a(i).*x0(i+1).^2+e(i)-c(i); W4(i)=a(i).*x0(i+1).^3-b(i).*x0(i).^3+c(i).*x0(i+1)-e(i).*x0(i); Si(t)=W1(i).*t^3+W2(i).*t^2+W3(i).*t+W4(i)%每个小区间的三次样条差值 函数表达式 end m=length(xx);T=zeros(m,1); fork=1:m forj=1:n if((xx(k)>=x0(j))&(xx(k) T(k)=W1(j).*(xx(k).^3)+W2(j).*(xx(k).^2)+W3(j).*xx(k)+W4(j); end end end T end %计算实习，第一种边界条件 clear; x=-1:0.2:1;%输入节点横坐标 y=ones(1,11)./(ones(1,11)+25*x.^2)%计算节点纵坐标 t=sym('t'); f=1/(1+25*t^2);%定义函数 f2=diff(f,2)%函数的二阶导数式 f21=vpa(subs(f2,'t',-1))%计算左端点的二阶导数值 . ;. f22=vpa(subs(f2,'t',1))%计算右端点的二阶导数值 xx=-1:0.1:1; T=mspline1(x,y,f21,f22,xx); T=T' ezplot(f,[-11]);%画出函数f的曲线 holdon plot(x,y,':',xx,T,'--');%根据函数计算和插值计算的结果画出的曲线 %计算实习，第二种边界条件 clear; x=-1:0.2:1;%输入节点横坐标 y=ones(1,11)./(ones(1,11)+25*x.^2)%计算节点纵坐标 t=sym('t'); f=1/(1+25*t^2);%定义函数 f1=diff(f,1)%函数的一阶导数式 f11=vpa(subs(f1,'t',-1))%计算左端点的一阶导数值 f12=vpa(subs(f1,'t',1))%计算右端点的一阶导数值 xx=-1:0.1:1; T=mspline2(x,y,f11,f12,xx); T=T' ezplot(f,[-11]);%画出函数f的曲线 holdon plot(x,y,':',xx,T,'--');%根据函数计算和插值计算的结果画出的曲线 4.计算结果 第一种边界条件： x-1-0.8-0.6-0.4-0.200.20.40.60.81 函数计算 值 0.038 5 0.058 8 0.10.20.510.50.20.1 0.058 8 0.038 5 插值计算 值 0.038 5 0.058 8 0.10.20.510.50.20.1 0.058 8 0 第二种边界条件： X-1-0.8-0.6-0.4-0.200.20.40.60.81 函数计算 值 0.03850.05880.10.20.510.50.20.10.05880.0385 . ;. 插值计算 值 0.03850.05880.10.20.510.50.20.10.05880