三次样条插值Word下载.docx
《三次样条插值Word下载.docx》由会员分享,可在线阅读,更多相关《三次样条插值Word下载.docx(11页珍藏版)》请在冰豆网上搜索。
![三次样条插值Word下载.docx](https://file1.bdocx.com/fileroot1/2022-10/24/a84f163d-d305-4a2b-a191-bb1ca3f8fcf0/a84f163d-d305-4a2b-a191-bb1ca3f8fcf01.gif)
(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