数值分析课程设计(三次样条插值)

《数值分析课程设计》

报告

专业：

学号：

学生姓名：

指导教师：

7.掌握三次样条插值函数的构造方法，体会三次样条插值函数对被逼近函数的近似。

三次样条插值函数边界条件由实际问题对三次样条插值在端点的状态要求给出。以第1 边界条件为例，

用节点处二阶导数表示三次样条插值函数，用追赶法求解相关方程组。通过Matlab 编制三次样条函数的通用程序，可直接显示各区间段三次样条函数体表达式，计算出已给点插值并显示各区间分段曲线图。

引言

分段低次样条插值虽然计算简单、稳定性好、收敛性有保证且易在电子计算机上实现，但只能保证各小段曲线在连接处的连续性，不能保证整件曲线的光滑性。利用样条插值，既可保持分段低次插值多项式，又可提高插值函数光滑性。故给出分段三次样条插值的构造过程算法步骤，利用Matlab软件编写三次样条插值函数通用程序，并通过数值算例证明程序的正确性。

三次样条函数的定义及特征

定义：设[a,b] 上有插值节点，a＝x1＜x2<…x n＝b，对应函数值为y1,y2,⋯y n。若函数S(x) 满足S(x j) ＝y j （j ＝1,2, ⋯,n ），S(x) 在[x j,x j+1] （j ＝1,2,⋯,n-1）上都是不高于三的多项式（为了与其对应j 从1 开始，在Matlab 中元素脚标从1 开始）。当S(x) 在[a,b] 具有二阶连续导数。则称S(x) 为三次样条插值函数。要求S(x) 只需在每个子区间[x j,x j+1] 上确定1 个三次多项式，设为：

Sj(x)＝ajx3＋bjx2＋cjx＋dj, (j＝1,2,⋯,n-1) (1)

其中a j,b j,c j,d j 待定，并要使它满足：

S(x j)＝y j, S(x j-0)=S(x j＋0), (j=2,⋯,n-1) (2)

S'(x j-0)=S'(x j＋0), S"(x j-0)=S"(x j＋0), (j=2,⋯,n-1) (3)

式(2)、(3)共给出n＋3(n-2)＝4n-6 个条件，

需要待定4(n-1) 个系数，因此要唯一确定三次插值函数，还要附加2 个边界条件。通常由实际问题对三次样条插值在端点的状态要求给出。常用边界的条件有以下3 类。

第1 类边界条件：给定端点处的一阶导数值，S'(x1)＝y1'，S'(x n)＝y n'。

第 2 类边界条件：给定端点处的二阶导数值，S"(x1)＝y1"，S"(x n)＝y n"。特殊情况y1"＝y n"＝0，称为自然边界条件。

第3 类边界条件是周期性条件，如果y＝f(x)是以b-a 为周期的函数，于是S(x) 在端点处满足条件S'(x1＋0)＝S'(x n-0)，S"(x＋0)＝S"(x n-0)。

下以第1 边界条件为例，利用节点处二阶导数来表示三次样条插值函数，给出具体的推导过程。

2 三次样条插值函数的推导过程

注意到S(x) 在[x j, x j+1]（j=1,2,⋯,n-1）上是三次多项式，于是S"(x)在[x j, x j+1] 上是一次多项式，如果S"(x) 在[x j,x j+1]（j=1,2,⋯,n-1）两端点上的值已知，设S"(x j)=M j，S"(x j+1)=M j+1，其中h j ＝x j+1-x j，对S"(x) 进行两次积分，则得到1 个具有2个任意常数A j，B j 的S(x) 表达式。对S"(x) 求两次积分

4 三次样条插值函数的Matlab 程序设计在Matlab[3]环境下根据上述算法步骤进行编程，源程序如下：

function []=spline3(X,Y,dY,x0,m)

N=size(X,2);

s0=dY(1); sN=dY(2);

interval=0.025;

disp('x0为插值点')

x0

h=zeros(1,N-1);

for i=1:N-1 h(1,i)=X(i+1)-X(i); end

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)); endmu=zeros(1,N-1); md=zeros(1,N-1);

md(1,N-1)=1; mu(1,1)=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

pp(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)&x0

y0=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

clear val fval

ans=sym(coeff,'d');