三次样条插值.docx
- 文档编号:13858163
- 上传时间:2023-06-18
- 格式:DOCX
- 页数:13
- 大小:136KB
三次样条插值.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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 三次 样条插值