三次样条插值方法的应用

CENTRAL SOUTH UNIVERSITY 数值分析实验报告

三次样条插值方法的应用

一、问题背景

分段低次插值函数往往具有很好的收敛性，计算过程简单，稳定性好，并且易于在在电子计算机上实现，但其光滑性较差，对于像高速飞机的机翼形线船体放样等型值线往往要求具有二阶光滑度，即有二阶连续导数，早期工程师制图时，把富有弹性的细长木条（即所谓的样条）用压铁固定在样点上，在其他地方让他自由弯曲，然后沿木条画下曲线，称为样条曲线。样条曲线实际上是由分段三次曲线并接而成，在连接点即样点上要求二阶导数连续，从数学上加以概括就得到数学样条这一概念。下面我们讨论最常用的三次样条函数及其应用。

二、数学模型

样条函数可以给出光滑的插值曲线（面），因此在数值逼近、常微分方程和偏微分方程的数值解及科学和工程的计算中起着重要的作用。

设区间[]b ,a 上给定有关划分b x x n =<<<= 10x a ，S 为[]b ,a 上满足下面条件的函数。

● )(b a C S ,2∈；

● S 在每个子区间[]1,+i i x x 上是三次多项式。

则称S 为关于划分的三次样条函数。常用的三次样条函数的边界条件有三种类型：

● Ⅰ型 ()()n n n f x S f x S ''0'',==。

● Ⅱ型 ()()n n n f x S f x S ''''0'''',==，其特殊情况为()()0''''==n n x S x S 。 ● Ⅲ型 ()() 3,2,1,0,0==j x S x S n j j ，此条件称为周期样条函数。

鉴于Ⅱ型三次样条插值函数在实际应用中的重要地位，在此主要对它进行详细介绍。

三、算法及流程

按照传统的编程方法，可将公式直接转换为MATLAB 可是别的语言即可；另一种是运用矩阵运算，发挥MATLAB 在矩阵运算上的优势。两种方法都可以方便地得到结果。方法二更直观，但计算系数时要特别注意。这里计算的是方法一的程序，采用的是Ⅱ型边界条件，取名为spline2.m 。

Matlab 代码如下：

function s=spline2(x0,y0,y21,y2n,x)

%s=spline2(x0,y0,y21,y2n,x)

%x0,y0 are existed points,x are insert points,y21,y2n are the second %dirivitive numbers given.

n=length(x0);

km=length(x);

a(1)=-0.5;

b(1)=3*(y0(2)-y0(1))/(2*(x0(2)-x0(1)));

for j=1:(n-1)

h(j)=x0(j+1)-x0(j);

end

for j=2:(n-1)

alpha(j)=h(j-1)/(h(j-1)+h(j));

beta(j)=3*((1-alpha(j))*(y0(j)-y0(j-1))/h(j-1)+alpha(j)*(y0(j+1)-y0(j ))/h(j));

a(j)=-alpha(j)/(2+(1-alpha(j))*a(j-1));

b(j)=(beta(j)-(1-alpha(j))*b(j-1))/(2+(1-alpha(j))*a(j-1)); end

m(n)=(3*(y0(n)-y0(n-1))/h(n-1)+y2n*h(n-1)/2-b(n-1))/(2+a(n-1)); for j=(n-1):-1:1

m(j)=a(j)*m(j+1)+b(j);

end

for k=1:km

for j=1:(n-1)

if ((x(k)>=x0(j))&(x(k)

l(k)=j;

end

end

end

for k=1:km

sum=(3*(x0(l(k)+1)-x(k))^2/h(l(k))^2-2*(x0(l(k)+1)-x(k))^3/h(l(k))^3)*y0(l(k));

sum=sum+(3*(x(k)-x0(l(k)))^2/h(l(k))^2-2*(x(k)-x0(l(k)))^3/h(l(k))^3)*y0(l(k)+1);

sum=sum+h(l(k))*((x0(l(k)+1)-x(k))^2/h(l(k))^2-(x0(l(k)+1)-x(k))^3/h(l(k))^3)*m(l(k));

s(k)=sum-h(l(k))*((x(k)-x0(l(k)))^2/h(l(k))^2-(x(k)-x0(l(k)))^3/h(l(k ))^3)*m(l(k)+1);

end

四、计算结果及分析

给定如下数值表，试求三次样条插值函数满足边界条件()()030''7.28''==s S

x28.7 28 29 30 ()x f 4.1 4.3 4.1 3.0

在MATLAB命令窗口中输入：

x=[28.7 28 29 30];

y=[4.1 4.3 4.1 3.0];

x0=[28.7:0.15:30];

y=spline2(x,y,0,0,x0);

plot(x0,y)

得到的图形如下图所示：