spline插值.docx

例5.6.1 给定以下数据, 求出三次样条函数，并计算函数分别在-0.15，-0.05，0.05，0.18，0.25处的近似值，并作图。
x
0.1
0.2
0.15
0
-0.2
0.3
y
0.95
0.84
0.86
1.06
0.50
0.72
解：编程如下:
clear
x=[0.1,0.2,0.15,0,-0.2,0.3];y=[0.95,0.84,0.86,1.06,1.50,0.72];
pp=spline(x,y);
pp.coefs
xx=[-0.15,-0.05,0.05,0.18,0.25];
yy=ppval(pp,xx) %or:yy=spline(x,y,xx)
fnplt(pp,'k')
hold on
plot(x,y,'o')
hold on
plot(xx,yy,'r*')
运行结果:
ans =
-36.3850 21.8592 -5.1164 1.5000
-36.3850 0.0282 -0.7390 1.0600
227.6995 -10.8873 -1.8249 0.9500
-143.0047 23.2676 -1.2059 0.8600
-143.0047 1.8169 0.0484 0.8400
yy =
1.2943 1.1016 1.0186 0.8409 0.8291

写出三次样条函数为:

练习: 给定函数，，节点，用三次样条插值求及在某些点处的近似值.

MATLAB实现:
yi=interp1(x,y,xi) %根据数据(x,y)给出分段线性插值函数的xi的值yi.
yi=interp1(x,y,xi,’spline’) %使用三次样条插值.
yi=interp1(x,y,xi,’cubic’) %使用分段三次插值.
Yi=spline(x,y,xi) %同yi=interp1(x,y,xi,’spline’)
pp=spline(x,y) %返回三次样条插值的分段多项式(pp)形式结构(使用“非扭结”端点条件).
MATLAB中的spline函数可以对数据点进行三次插值，默认的边界条件强制使得插值第一段的三次项系数与第二段的三次项系数相同，最后一段与倒数第二段的多项式系数相同，这就所谓的非扭结（not-a-knot）条件.
1一维插值函数
 Interp1（）
命令格式：
yi=interp1（x,y,xi,’method’）
x为插值节点构成的向量，y为插值节点函数值构成的向量，yi是被插值点xi的插值结果，‘method‘是采用的插值方法，缺省时表示分线段性插值，’nearest‘为最邻近插值；’linear‘为分线段性插值;’spline’为三次样条插值；’pchip’为分段Hermite插值；’cubic’为分段Hermite插值
例子：画出y=sin(x)在区间[0 10]的曲线，并在曲线上插值节点xk=k，k=0,1…10及函数值，画出分段线性插值折线图
x=0:10;
y=sin(x);
xi=0:0.25:10;
yi1=interp1(x,y,xi,'nearest');
yi2=interp1(x,y,xi,'linear');
yi3=interp1(x,y,xi'spline');
yi4=interp1(x,y,xi,'pchip');
yi5=interp1(x,y,xi,'cubic');
subplot(1,5,1)
plot(x,y,'o',xi,yi1,'k--',xi,sin(xi),'k:');
title('\bfNearest');
subplot(1,5,2)
plot(x,y,'o',xi,yi2,'k--',xi,sin(xi),'k:');
title('\bfLinear');
subplot(1,5,3)
plot(x,y,'o',xi,yi3,'k--',xi,sin(xi),'k:');
title('\bfSpline');
subplot(1,5,4)
plot(x,y,'o',xi,yi4,'k--',xi,sin(xi),'k:');
title('\bfPchip');
subplot(1,5,1)
plot(x,y,'o',xi,yi5,'k--',xi,sin(xi),'k:');
title('\bfCubic');

spline()为三次样条函数
命令格式1：yi=spline(x,y,xi),意义等同于yi=interp1(x,y,xi,'spline')
命令格式2：pp=spline(x,y) ，输出三次样条函数分段表示的结构
pchip()命令格式与spline()完全相同

csape()为可输入边界条件的三次样条函数
命令格式：pp=csape(x,y,conds,valconds),x为插值节点构成的向量，y为插值节点函数值构成的向量；conds为边界类型，缺省为非扭结边界条件；valconds表示边界值。
边界类型：‘complete‘为给定边界条件的一阶导数；’not-a-knot‘为非扭结；’periodic‘为周期边界条件；’second’为给定边界条件的二阶导数;’variational’为自由边界
例子：
（1）已知y=f(x)的函数表及端点条件S’’(x1)=S’’(x4)=0
X     1     2     4     5
F(x)       1     3     4     2
求三次样条插值函数S（x）,并计算f（3），f(4.5)的近似值
clear;
clc;
x=[1 2 4 5];
y=[1 3 4 2];
s=csape(x,y,'variation')
value=fnval(s,[3 4.5])
再输入s.coefs可以得到三次样条插值分段表示的系数
注释：V = FNVAL(F,X)  or FNVAL(X,F)  provides the value atthe points in  X  of the function described by  F .
（1）已知函数y=1/(25x^2+1)在[0 1]上的值如下表
X              0     0.25            0.5        0.75            1
Y              1  0.3903       0.1379       0.0664       0.0385
求三次样条插值函数S(x)，使满足S’(0)=0,S’(1)=-0.074
x=[0:0.25:1];
y=1./(1+25*x.^2);
s=csape(x,y,'complete',[0 -0.074])
fnplt(s,'r')
s.coefs
题目如下：
清华大学出版社的《数值分析（第5版）》
P49,20题。
x=[0.25 0.3 0.39 0.45 0.53];
y=[ 0.5 0.5477 0.6245 0.6708 0.7280 ]
pp=csape(x,y,'second',[0,0.0]);
disp(pp.coefs);

其中COEFS的含义是在Xi-Xi+1区间上的多项式是，例如COEFS数组第一行的意思是在X=0.25到X=0.3的区间上时表达式是-6.2652*（X-0.25)^3+0.9697*(X-0.25)^1+0.5;

   -6.2652    0.0000    0.9697    0.5000
    1.8813   -0.9398    0.9227    0.5477
   -0.4600   -0.4318    0.7992    0.6245
    2.1442   -0.5146    0.7424    0.6708
csape，是计算在各种边界条件下的三次样条插值。
pp = csape(x,y,conds)
其中conds主要有以下的选项variational（自然边界条件，首末点二阶导数均为0），second（指定首末点的二阶导数），periodic（周期性边界条件，首末点的0~2阶导数相等），complete（给定导数情况，默认）

function pp = csape(x,y,conds,valconds)
%pp=csape(x,y,'变界类型','边界值'),生成各种边界条件的三次样条插值. 其中,(x,y)为数据向量
%边界类型可为:'complete',给定边界一阶导数.
%             'not-a-knot',非扭结条件,不用给边界值.
%             'periodic',周期性边界条件,不用给边界值.
%             'second',给定边界二阶导数.
%             'variational',自然样条(边界二阶导数为0)
% .
%例 考虑数据
%          x | 1 2 4 5
%         ---|-------------
%          y | 1 3 4 2
%边界条件S''(1)=2.5,S''(5)=-3,
%     x=[1 2 4 5];y=[1 3 4 2];
%     pp=csape(x,y,'second',[2.5,-3]);pp.coefs
%      xi=1:0.1:5;yi=ppval(pp,xi);
%      plot(x,y,'o',xi,yi);
matlab中的插值问题转自mop)
matlab中的插值问题
今天写数值分析作业，仔细的研究了一下matlab中的插值问题，现总结如下：
interp1 ， 是一维数据的插值函数，基本使用方法如下
yi = interp1(x,Y,xi,method)
注意，这里常用的method有linear,spline。其中linear（线性）是默认的方法。spline应该和使用spline是一样的。
spline，三次样条插值。注意它默认使用的是not-a-knot边界条件，也就是第一个点的三次导数和第二点的三次导数一样；最后一个点的三次导数和倒数第一个点一样。当y=[df1,y,df2]时，表示第一点和第二个点的一阶导数分别为df1,df2。
x = -4:4;
y = [0 .15 1.12 2.36 2.36 1.46 .49 .06 0];
cs = spline(x,[0 y 0]);
xx = linspace(-4,4,101);
plot(x,y,'o',xx,ppval(cs,xx),'-');
x = pi*[0:.5:2]; 
y = [0 1 0 -1 0 1 0; 
1 0 1 0 -1 0 1];
pp = spline(x,y);
yy = ppval(pp, linspace(0,2*pi,101));
(以上例子参考matlab的help)

上例中pp是一种多项式的表达方式，通过ppval就能求解出相应的值。
csape，是计算在各种边界条件下的三次样条插值。
pp = csape(x,y,conds)
其中conds主要有以下的选项variational（自然边界条件，首末点二阶导数均为0），second（指定首末点的二阶导数），periodic（周期性边界条件，首末点的0~2阶导数相等），complete（给定导数情况，默认）

function pp = csape(x,y,conds,valconds)
%pp=csape(x,y,'变界类型','边界值'),生成各种边界条件的三次样条插值. 其中,(x,y)为数据向量
%边界类型可为:'complete',给定边界一阶导数.
%             'not-a-knot',非扭结条件,不用给边界值.
%             'periodic',周期性边界条件,不用给边界值.
%             'second',给定边界二阶导数. 
%             'variational',自然样条(边界二阶导数为0)
% .
%例 考虑数据
%          x | 1 2 4 5
%         ---|-------------
%          y | 1 3 4 2
%边界条件S''(1)=2.5,S''(5)=-3,
%     x=[1 2 4 5];y=[1 3 4 2];
%     pp=csape(x,y,'second',[2.5,-3]);pp.coefs
%      xi=1:0.1:5;yi=ppval(pp,xi);
%      plot(x,y,'o',xi,yi);
pp0 = csape(x,[1,zeros(1,length(y)),0],[1,0]);
pp = csape( x, [1 sin(x) 0], [1 2] ) ％左边的点一阶导数为1，右边的点二阶导数为0
splinetool是一个图形化的插值工具
lagrange插值，由于lagrange插值可能不收敛，所以工程中很少有人用这种插值。matlab中没有专门的lagrange插值函数。但我们可以自己编一个，如下：
%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

6