第六章 插值与拟合

y(x)
h [( x
i i 1
n
i
x )( 2 a i y i y i ) y i ]
其中
n
2.MATLAB实现
function y=hermite(x0,y0,y1,x) %hermite interpolation n=length(x0);m=length(x); for k=1:m yy=0.0; for i=1:n h=1.0; a=0.0; for j=1:n if j~=i h=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2; a=1/(x0(i)-x0(j))+a; end end yy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i)); end y(k)=yy; end
1.方法介绍 在科学实验的统计方法研究中，往往要从一组实验 数据 ( x i , y i )中寻找出自变量x和因变量y之间的函数关 系y=f(x)。由于观测数据往往不够准确，因此并不要 ) 求y=f(x)经过所有的点 ( x i , y i ，而只要求在给定点 x i 上误差 i f ( x i ) y i 。按照某种标准达到最小，通常 采其欧氏范数 作为误差量度的标准。这就是所 谓的最小二乘法。
1.方法介绍 不少实际问题中不但要求在节点上函数值相等，而 且要求导数值也相等，甚至要求高阶导数值也相等， 满足这一要求的插值多项式就是Hermite插值多项 式。下面只讨论函数值与一阶导数值个数相等且已 知的情况。 已知n个插值节点 x 1 , x 2 , , x n 及其对应的函数 值 y 1 , y 2 , , y n 。和一阶导数值 y 1 , y 2 , , y n 。 y 即： ( x i ) y i , y ( x i ) y i , i 1, 2 , , n 。则计算插值 区域内任意x的函数值y的Hermite插值公式为:
y(x)

n
k 1
yk
n
x xj xk x j
j 1 jk
2.MATLAB实现
%lagrange interpolation 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
所有的插值方法要求x0是单调的。x0也可能并非连续 等距的。当x0为等距且连续时可以用快速插值法。 使用这些方法的格式分别为'*1inear'、'*cubic'、 '*nearest'或'*spline'。x0为非连续的快速线性插值的 情况interplq。
【例6.1.3】正弦曲线的插值示例。
>>x0=0:1:10; y0=sin(x0); x1=0:0.25:10; y1=interp1(x0,y0,x1); plot(x0,y0,'o',x1,y1)
0.6
同类的函数 还有 Interplq interpfl spline interp2 interpn。
0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1
0
1
2
3
4
5
6
7
8
9
10
>> x0=[-5:1:5]; >> y0=1./(1+x0.^2); >> x=[-5:0.1:5]; >> y=interp1(x0,y0,x); >> y1=1./(1+x.^2); >> plot(x,y,'--r',x,y1, '-b') >> legend('插值曲线', '原曲线')
x0=[-5:1:5]; y0=1./(1+x0.^2); x=[-5:0.1:5]; y=interp1(x0,y0,x,'spline'); y1=1./(1+x.^2); plot(x,y,'--r',x,y1, '-b') legend('插值曲线', '原曲线')
6.1.3 Hermite插值
3.习题举例
【例6.2.1】设 下表所示的数据。
y a b x cx ，用最小二乘法拟合如
2
x 0.5 1.0 1.5 2.0 2.5 3.0 y 1.75 2.45 3.81 4.80 8.00 8.60 此例用polyfit功能函数进行拟合。 在MATLAB 命令窗中输入
x0=[0.5 1.0 1.5 2.0 2.5 3.0]; y0=[1.75 2.45 3.81 4.80 8.00 8.60]; a=polyfit(x0,y0,2)； x1=[0.5:0.05:3.0]; y1=a(1)*x1.^2+a(2)*x1+a(3); plot(x0,y0,'*') hold on plot(x1,y1,'-r')
10 9 8 7 6 5 4 3 2 1 0.5
1
1.5
2
2.5
3
【例6.2.2】用最小二乘法求一个形如 y a bx 2 的经 验公式，使它与下表所示的数据相拟合。
xi yi 19 25 31 38 44
19.0 32.3 49.0 73.3 98.8
下面用另一种方法来求解此拟合问题。用求解矩 阵的方法来解，把a,b看成是未知量。已知xi,yi求 解一超定方程。在MATLAB中实现为: x0=[19 25 31 38 44]; y0=[19.0 32.3 49.0 73.3 98.8]; x1=x0.^2; x1=[ones(5,1),x1']; ab=x1\y0'; x=[19:0.2:44]; y=ab(1)+ab(2)*x.^2; plot(x0,y0,'o') hold on plot(x,y,'-r')
100
90
80
70
60
50
40
30
20
10 15
20
25
30
35
40
45
除了polyfit命令实现拟合外，matlab还可以通过 basic—fitting完成拟合的功能。在basic—fitting 中首先画出历史数据的散点图，分析其所呈现的规 律，再进行拟合。 例如： >> Year=[1625,1830,1930,1960,1974,1987,1999]; >> Population=[5,10,20,30,40,50,60]; >> plot(Year,Population,'X'); >> xlabel('年份') >> ylabel('人口数量（亿人）')
1 0.8 0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
0
1
2
3
4
5
6
7
8
9
10
x0=0:1:10; y0=sin(x0); x1=0:0.005:10; y1=interp1(x0,y0,x1,'nearest'); figure(2) 1 plot(x0,y0,'o',x1,y1) 0.8
>>x0=[-5:1:5]; y0=1./(1+x0.^2); x=[-5:0.1:5]; y=lagrange(x0,y0,x); y1=1./(1+x.^2); plot(x,y,'--r',x,y1, '-b') legend('插值曲线', '原曲线')
插值曲线 原曲线
为了解决Runge问题，引入分段线性插值。
hi

(
x xj xi x j
)
2
j 1 ji
ai

n
1 xi x j
j 1 i j
3.习题举例 x0=[0.3 0.32 0.35]; y0=[0.29552 0.31457 0.34290]; y1=[0.95534 0.94924 0.93937]; x=[0.3:0.005:0.35]; y=hermite(x0,y0,y1,x); plot(x,y) y=hermite(x0,y0,y1,0.34) y2=sin(x); hold on plot (x,y2,'--r')
来拟合这些数据点，解决步骤如下：
来拟合这些数据点，解决步骤如下：
1.分别输入数据矩阵和函数值向量。
x 11 x 21 ... xm1 x 12 x 22 ... xm 2 ... ... ... ... x1n x2n ... xmn Y0 y1 y2 ... ym
yl=interp1(y0,x) 此格式默认x0=l:n，n为向量y0的元素 个数值，或等于矩阵y0的size(y,1)。 yl=interp1(x0,y0,x,'Method')
method用来指定插值的算法。默认为线性算法,其值可 以是如下的字符串。 · nearest 线性最近项插值。 · linear 线性插值。 · spline 三次样条插值。 · cubic 三次插值