三次样条拟合范例

1设计目的、要求

对龙格函数2

2511

)(x

x f +=

在区间[-1,1]上取10=n 的等距节点，分别作多项式插值、三次样条插值和三次曲线拟合，画出)(x f 及各逼近函数的图形，比较各结果。

2设计原理

（1） 多项式插值：利用拉格朗日多项式插值的方法，其主要原理是拉格朗日多项

式，即：

01,,...,n x x x 表示待插值函数的1n +个节点，

()()n

n j k k j j k L x y l x y ===∑，其中0,1,...,j n =；

011011()...()()...()

()()...()...()...()

k k n k k k k k k k n x x x x x x x x l x x x x x x x x x -+-+----=

----

（2） 三次样条插值：三次样条插值有三种方法，在本例中，我们选择第一边界条

件下的样条插值，即两端一阶导数已知的插值方法：

00'()'S x f = '()'n n S x f =

（3）三次曲线拟合：本题中采用最小二乘法的三次多项式拟合。最小二乘拟合是

利用已知的数据得出一条直线或者曲线，使之在坐标系上与已知数据之间的距离的平方和最小。在本题中，n= 10，故有11个点，以这11个点的x 和

y 值为已知数据，进行三次多项式拟合，设该多项式为

23432xi i i i p a a x a x ax =+++，该拟合曲线只需2[]xi i p y -∑的值最小即可。

3采用软件、设备

计算机、matlab 软件

4设计内容

1、多项式插值：

在区间[]1,1-上取10=n 的等距节点，带入拉格朗日插值多项式中，求出各个节点的插值，并利用matlab 软件建立m 函数，画出其图形。

在matlab 中建立一个lagrange.m 文件，里面代码如下： %lagrange 函数

function y=lagrange(x0,y0,x) n=length(x0);m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j~=k

p=p*(z-x0(j))/(x0(k)-x0(j)); end end

s=p*y0(k)+s; end y(i)=s; end

建立一个polynomial.m 文件，用于多项式插值的实现，代码如下： %lagrange 插值 x=[-1:0.2:1];

y=1./(1+25*x.^2); x0=[-1:0.02:1];

y0=lagrange(x,y,x0); y1=1./(1+25*x0.^2); plot(x0,y0,'--r') %插值曲线 hold on %原曲线

plot(x0,y1,'-b')

运行duoxiangshi.m 文件，得到如下图形：

2、三次样条插值：

所谓三次样条插值多项式()n S x 是一种分段函数，它在节点

i x 011()n n a x x x x b -=<<⋅⋅⋅<<=分成的每个小区间1[,]i i x x -上是3次多项式，其在

此区间上的表达式如下：

2233

1111111()[()()]()()666[,]1,2,,.

i i i i i i i i i i i i i i i

i i h x x h x x S x x x M x x M y M y M h h h x x x i n --------=-+-+-+-∈=⋅⋅⋅，

， 因此，只要确定了i M 的值，就确定了整个表达式，i M 的计算方法如下： 令：

11111111116()6(,,)i i i i i i i i i i i i i i

i i i i i i i h h h h h h y y y y d f x x x h h h h μλμ++++--+++⎧

===-⎪++⎪⎨

--⎪=-=⎪+⎩，

则i M 满足如下1n -个方程：

1121,2,,1

i i i i i i M M M d i n μλ-+++==⋅⋅⋅-，

对于第一种边界条件下有

⎪⎪⎩⎪

⎪⎨

⎧

-=+-=+---00011

0111

)

'],([62]),['(62h f x x M M h x x f f M M n n n n n n

如果令10010001

6([,]')6('[,])

1,,1,,n n n n n n f x x f f f x x d d h h λμ----==

==那么解就可以为

⎪⎪

⎪

⎪⎪⎪⎭

⎫ ⎝⎛=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎪⎭⎫ ⎝

⎛----n n n n n n n d d d d M M M M 11011011

1

10

2222

μλμλμλ

求函数的二阶导数： >> syms x

>> f=sym(1/(1+25*x^2)) f =

1/(1+25*x^2) >> diff(f) ans =

-(50*x)/(25*x^2 + 1)^2

将函数的两个端点，代入上面的式子中：

f’(-1)=0.0740

f’(1)=-0.0740

求出从-1到1的n=10的等距节点，对应的x，y值

对应m文件代码如下：

for x=-1:0.2:1

y=1/(1+25*x^2)

end

y =

得出

x=-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

y=0.0385 0.0588 0.1 0.2 0.5 1 0.5 0.2 0.1 0.0588 0.0385

编写m文件Three_Spline.m

x=linspace(-1,1,11);

y=1./(1+25*x.^2);

[m,p]=scyt1(x,y,0.0740,-0.0740);

hold on

x0=-1:0.01:1;

y0=1./(1+25*x0.^2);

plot(x0,y0,'--b')

得到如下图像：