三次样条插值.docx
《三次样条插值.docx》由会员分享，可在线阅读，更多相关《三次样条插值.docx（13页珍藏版）》请在冰点文库上搜索。
三次样条插值
三次样条插值
三次样条插值
1.算法原理
由于在许多实际问题中，要求函数的二阶导数连续，人们便提出了三次样条插值函数，三次样条插值函数是由分段三次函数拼接而成的，在连接点处二阶导数连续。
设S（x）在节点
处的二阶导数
，其中
为待定参数。
由S（x）是分段三次多项式可知，
是分段线性函数，
在子区间
上可以表示为
其中
，对S（x）两端积分两次得
其中
和
为积分常数。
由插值条件
得
由此解得
代入得：
三种边界条件的三弯矩方程：
（1）第一种边界条件：
已知
。
取
，这时方程组减少了两个未知量，变成只含n-1个未知量
的n-1个方程的方程组
，用矩阵表示为
可用追赶法求解出
后，即得三条样条插值函数。
（2）第二种边界条件，已知
，记
，则有
，得
，即
，其中
，得到方程组
，用矩阵表示为
，该方程组的系数矩阵是严格三对角占优矩阵，可用追赶法求解。
（3）第三种边界条件：
周期型边界条件。
已知
是以
为周期的周期函数，则由周期性可知，
，这时将点
看成是内节点，则有
，也即
，其中
，方程组第1个方程为：
，所以方程组为
，用矩阵表示为
，显然系数矩阵为严格对角占优矩阵，可用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.2
0
0.2
0.4
0.6
0.8
1
函数计算值
0.0385
0.0588
0.1
0.2
0.5
1
0.5
0.2
0.1
0.0588
0.0385
插值计算值
0.0385
0.0588
0.1
0.2
0.5
1
0.5
0.2
0.1
0.0588
0
第二种边界条件：
X
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
函数计算值
0.0385
0.0588
0.1
0.2
0.5
1
0.5
0.2
0.1
0.0588
0.0385
插值计算值
0.0385
0.0588
0.1
0.2
0.5
1
0.5
0.2
0.1
0.0588
0
升级会员 