三次样条插值Word下载.docx
《三次样条插值Word下载.docx》由会员分享，可在线阅读，更多相关《三次样条插值Word下载.docx（11页珍藏版）》请在冰豆网上搜索。
（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:
y（i）=d（i）-l（i）*y（i-1）;
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）;
x
functionT=mspline1（x0,y0,f21,f22,xx）
%三次样条插值函数第一种边界条件
%x0、y0分别为节点的横坐标和纵坐标；
%f21为左端点的二阶导数值,f22为右端点的二阶导数值；
xx为由插值点组成的向量
n=length（x0）-1;
%计算小区间数
fori=1:
h（i）=x0（i+1）-x0（i）;
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））;
A=zeros（n-1）;
n-2
A（i+1,i）=mu（i+1）;
%次下对角线
A（i,i+1）=lamda（i）;
%次上对角线
A（i,i）=2;
%主对角线
A（n-1,n-1）=2;
dd=zeros（n-1,1）;
%右端列向量
dd（i）=d（i）;
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）;
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）%每个小区间的三次样条差值函数表达式
m=length（xx）;
T=zeros（m,1）;
fork=1:
m
forj=1:
if（（xx（k）>
=x0（j））&
（xx（k）<
x0（j+1）））
T（k）=W1（j）.*（xx（k）.^3）+W2（j）.*（xx（k）.^2）+W3（j）.*xx（k）+W4（j）;
end
T
End
functionT=mspline2（x0,y0,f11,f12,xx）
%三次样条插值函数第二种边界条件
%f11为左端点的二阶导数值,f12为右端点的二阶导数值；
A=zeros（n+1）;
%系数矩阵
dd=zeros（n+1,1）;
A（k,k）=2;
%主对角线元素
A（k,k-1）=mu（k-1）;
%次下对角线元素
A（k,k+1）=lamda（k-1）;
%次上对角线元素
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）;
dd（k）=d（k-1）;
dd
M
%计算实习，第一种边界条件
clear;
x=-1:
0.2:
1;
%输入节点横坐标
y=ones（1,11）./（ones（1,11）+25*x.^2）%计算节点纵坐标
）;
f=1/（1+25*t^2）;
%定义函数
f2=diff（f,2）%函数的二阶导数式
f21=vpa（subs（f2,'
-1））%计算左端点的二阶导数值
f22=vpa（subs（f2,'
1））%计算右端点的二阶导数值
xx=-1:
0.1:
T=mspline1（x,y,f21,f22,xx）;
T=T'
ezplot（f,[-11]）;
%画出函数f的曲线
holdon
plot（x,y,'
:
'
xx,T,'
--'
%根据函数计算和插值计算的结果画出的曲线
%计算实习，第二种边界条件
f1=diff（f,1）%函数的一阶导数式
f11=vpa（subs（f1,'
-1））%计算左端点的一阶导数值
f12=vpa（subs（f1,'
1））%计算右端点的一阶导数值
T=mspline2（x,y,f11,f12,xx）;
4.计算结果
第一种边界条件：
-1
-0.8
-0.6
-0.4
-0.2
0.2
0.4
0.6
0.8
函数计算值
0.0385
0.0588
0.1
0.5
插值计算值
第二种边界条件：
X
升级会员 