三次样条插值的求解

三次样条插值的求解

摘要：分段低次插值虽然解决了高次插值的振荡现象和数值不稳定现象，使得插值多项式具有一致收敛性，保证了插值函数整体的连续性，但在函数插值节点处不能很好地保证光滑性要求，这在某些要求光滑性的工程应用中是不能接受的。如飞机的机翼一般要求使用流线形设计，以减少空气阻力，还有船体放样等的型值线，往往要求有二阶光滑度（即有二阶连续导数）。因此，在分段插值的基础上，引进了一种新的插值方法，在保证原方法的收敛性和稳定性的同时，又使得函数具有较高的光滑性的样条插值。

关键字:三转角方程 三弯矩阵方程

0. 引言

1，三次样条函数

x,x定义1：若函数S(x)C[a,b],且在每个小区间上jj1上是三次多

2项式，其中ax0x1xnb 是给定节点，则称s(x)是节点

x0,x1,,xn上的三次样条函数。若节点xj上 给定函数值

yjf(xj)(j0,1,n) ,且

s(xj)yj (j0,1,n) （1.1）

成立，则称 s(x)为三次样条差值函数。

从定义知，要求出s(x)，在每个应小区间[xj,xj1] 上确定4个待定系数，共有 n个小区间，故应确定4n 个参数，根据s(x) 在[a,b] 上二阶导数连续，在节点xjj1,2,3,,n1处应满足连续性条件

s(xj0)s(xj0), s'(xj0)s'(xj0),

s''(xj0)s''(xj0) （1.2） 共有 3n-3个条件，再加上s(x)满足插值条件（1.1），共有4n-2个条件，因此还需要2个条件才能确定s(x)。通常可在区间[a,b] 端点ax0,bxn上各加一个条件（称边界条件），边界条件可根据实际的问题要求给定。常见的三种： (1) 已知两端的一节导数值，即

s(x0)f0''''s(xn)fn （1.3）

(2)两端的二阶导数已知，即

s(x0)f0s(xn)fn'''''''' （1.4）

特殊情况下的边界条件

s''(x0)s''(xn)0 （1.4）’ 称为自然边界条件

（3）当f(x)是以xnx0 为周期函数时，则要求s(x) 也是周期函数，这时边界条件应满足

而此时式中

2.三转角方程

及相应的边界条件函数s(x)的表达式，若满足假定的

在节点

, 这样确定的样条函数

称为周期函数。

处的值为,再由，是 ,则由分段Hermite插

值式

得 （2.1）

其中是插值基函数。由

得，（2.1）中

在整个区间

上连续。且 满足

。然而在

公式（2.1）中是未知的,它可利用式

及边界条件

来确定。为了求出，下面考虑是

这里 ,对

于是 s''(xj0)4hjmj2hjmj16h2j(yj1yj)

同理可得s''(x)在区间[xj1,xj]上的表达式

S(x)6(xj1xj2x)h3j1(yjyj1)6x2xj14xjh2j1mj16x4xj12xjh2j1mj

及 s''(xj0)由

4hj1mj2hj1mj16h2j1(yjyj!)

得

1hj1mj12(1hj11hj)mj1hjmj13(yj1yjh2jyjyj1h2j1)(j1,2,,n1)

用

1hj11hj除上式的两边,并加以整理,得

jmj12mjjmj1gj 其中：(j1,2,,n1)jghjhjhj1jhj1hjhj1yj1yjj（2.2）

)(j1,2,,n1)j3(jyjyj1hj1hj共个n1个方程,n1个未知量(2.2)式称为基本方程组

如果问题要求满足第一类(一阶)边界条件:

即 m0f0 mnfn

基本方程组(2.2)化为n-1阶方程组

3.算法步骤： ①

输 入 n 个插值结点，ax0x1xnb，对应的函数插值为

y1,y2.,yn,边界条件y1,y2''，待求插值点x0。

② ③

计算计

hjxjxj (j1,2,3,n1) 。

1hj1hj算

jhj1hj1hj ，

dj6(yj1yjhj6h1(yjyj1hj1'),(j2,,n1) 。

④ ⑤ ⑥ ⑦

计算1y2y1h1y1)， n6hn1(yn'ynyn1hn1)。

用追赶法求解方程组() 输出各区间的sj(x)表达式。

判断x0所在区间[xj,xj1]并计算插值y0，绘制各区间的三次样条插值曲线。

4.三次样条插值函数的Matlab程序设计

在Matlab环境下根据上述算法步骤进行编程，源程序如下： N=size( X,2);

function []=spline3(X,Y,dY,x0,m) s0=dY(1); sN=dY(2); inte rval=0.025; disp('x0为 插值 点') x0

for i=1:N-1 h(1,i) =X(i+1) -X(i); e nd h=ze ros(1,N-1) ;

d(1,1)=6*( (Y(1,2) -Y(1,1) )/h( 1,1)- s0)/h(1,1); d(N,1)=6*(sN-(Y(1,N)-Y(1,N-1))/h(1,N-1))/h(1,N-1); for i=2:N-1

d(i,1)=6*((Y(1,i+1)-Y(1,i))/h(1,i)-(Y(1,i)-Y(1,i-1))/h(1,i-1) ( h( 1,i)+h( 1,i-1)) ; end

mu=zer os(1,N-1); md=zer os( 1,N-1); md=zeros(1,N-1); for i=1:N-2

u=h(1,i+1)/(h(1,i)+h(1,i+1)); mu(1,i+1)=u;

md(1,i)=1-u; end

p(1,1)=2; q(1,1)=mu(1,1)/2; for i=2:N-1

p(1,i)=2-md(1,i-1)*q(1,i-1); q(1,i)=mu(1,i)/p(1,i); end p(1,N)=2-md(1,N-1)*q(1,N-1); y=zeros(1,N); y(1,1)=d(1)/2;

for i=2:N y(1,i)=(d(i)-md(1,i-1)*y(1,i-1))/p(1,i); end x=zeros(1,N); x(1,N)=y(1,N);

for i=N-1:-1:1 x(1,i)=y(1,i)-q(1,i)*x(1,i+1); end fprintf ('M为三对角方程的解\\n'); M=x; fprintf ('\\n'); syms t;

digits (m); for i=1:N-1

2pp(i)=M(i)*(X(i+1)-t)^3/(6*h(i))+M(i+1)*(t-X(i))^3 /(6*h(i))+(Y(i)-M(i)*h(i)^2/6)*(X(i+1)-t)/h(i)+ (Y(i+1)-M(i+1)*h(i)^2/6)*(t-X(i))/h(i);

pp(i)=simplify(pp(i)); coeff=sym2poly(pp(i)); if length(coeff)~=4

tt=coeff(1:3); coeff(1:4)=0; coeff(2:4)=tt; end if x0>X(i)&x0y0=coeff(1)*x0^3+coeff(2)*x0^2+coeff(3)*x0+coeff(4); end

val=X(i):interval:X(i+1); for k=1:length(val)

fval(k)=coeff(1)*val(k)^3+coeff(2)*val(k)^2

)/

+coeff(3)*val(k)+coeff(4); end

if mod(i,2)==1 plot(val,fval,'r+') else plot(val,fval,'b.') end hold on

ans=poly2sym(ans,'t')

fprintf('在区间[%f,%f]内\\n',X(i),X(i+1)); fprintf('三次样条函数S(%d)=',i); pretty(ans); end

fprintf ('x0所在区间为[%f,%f]\\n',X(L),X(L+1)); pretty(ans); end fprintf ('函数在插值点x0=%f的值为\\n',x0); y0

程序中：X,Y为输入结点，dY为两端点一阶导阵，x 为待求插值点，m为有效数字位数。应

y＝f(x)的函数值 x -1.5 0 1 2 y 0.125 -1 1 9

syms x

f1=1/(1+x^2); f2=1/(1+5*x^2); x=-5:5;

y1=subs(f1,x); y2=subs(f2,x); xx=-5:0.25:5;

yy=spline(x,y1,xx); yy2=spline(x,y2,xx); subplot(2,2,1);

plot(x,y1,'o'.xx.yy), grid

t=-5:0.25:5; yt=subs(f1,t); yt2=subs(f2,t); subplot(2,2,2); plot(t,yt). grid

subplot(2,2,3)

plot(x,y2,'o',xx,yy2), grid

subplot(2,2,4) plot(t,yt2), grid