[数学建模(六)]使用MATLAB实现插值

常用的插值：拉格朗日多项式插值、牛顿插值、分段线性插值、Hermite 插值和三次样条插值。

1.拉格朗日插值法

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=1:5;

>> y0=x0.^2;

>> x=2.5;

>> y=lagrange(x0,y0,x)

y =

6.2500

>> plot(x0,y0)

>> hold on

>> plot(x,y,'*')

2.牛顿插值

function yi=newton(x,y,xi)

n=length(x);

m=length(y);

if n~=m

    error('The lengths of X ang Y must be equal!');

    return;

end

Y=zeros(n);

Y(:,1)=y';

for k=1:n-1

    for i=1:n-k

        if abs(x(i+k)-x(i))<eps

            error('the DATA is error!');

            return;

        end

        Y(i,k+1)=(Y(i+1,k)-Y(i,k))/(x(i+k)-x(i));

    end

end

yi=0;

for i=1:n

    z=1;

    for k=1:i-1

        z=z*(xi-x(k));

    end

    yi=yi+Y(1,i)*z;

end

测试：

x0=1:5;

y0=x0.^2;

x=2.5;

y=newton(x0,y0,x)

plot(x0,y0)

hold on

plot(x,y,'*')

y =

6.2500

3.分段线性插值

y=interp1(x0,y0,x,'method')

method 指定插值的方法，默认为线性插值。其值可为：

'nearest' 最近项插值

'linear' 线性插值

'spline' 逐段3 次样条插值

'cubic' 保凹凸性3 次插值。

所有的插值方法要求 x0 是单调的。

测试：

x0=1:5;

y0=x0.^2;

x=2.5;

y=interp1(x0,y0,x,'linear')

plot(x0,y0)

hold on

plot(x,y,'*')

y =

6.5000

4. 埃尔米特(Hermite)插值

function y=hermite(x0,y0,y1,x); %y1是导数值

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

测试：

>> x0=1:5;

y0=x0.^2;

x=2.5;

y1=[1 1 1 1 1];

y=hermite(x0,y0,y1,x)

plot(x0,y0)

hold on

plot(x,y,'*')

y =

7.4228

5.样条插值

Matlab 中三次样条插值也有现成的函数：

(1)y=interp1(x0,y0,x,'spline');

(2)y=spline(x0,y0,x);

(3)pp=csape(x0,y0,conds);

y=ppval(pp,x)

说明：

csape 的返回值是pp 形式，要求插

值点的函数值，必须调用函数ppval。

pp=csape(x0,y0)：使用默认的边界条件，即Lagrange 边界条件。

pp=csape(x0,y0,conds)中的conds 指定插值的边界条件，其值可为：

'complete' 边界为一阶导数，即默认的边界条件

'not-a-knot' 非扭结条件

'periodic' 周期条件

'second' 边界为二阶导数，二阶导数的值[0, 0]。

'variational' 设置边界的二阶导数值为[0,0]。

测试：

x0=1:5;

y0=x0.^2;

x=2.5;

y1=interp1(x0,y0,x,'spline')

y2=spline(x0,y0,x)

pp=csape(x0,y0,'second')

y3=ppval(pp,x)

plot(x0,y0)

hold on

plot(x,y1,'*')

plot(x,y2,'*')

plot(x,y3,'*')

y1 =

6.2500

y2 =

6.2500

pp =

包含以下字段的 struct:

form: 'pp'

breaks: [1 2 3 4 5]

coefs: [4×4 double]

pieces: 4

order: 4

dim: 1

y3 =

6.2321

6.二维插值

一维插值：节点为一维变量，插值函数是一元函数（曲线）。

二维插值：节点是二维的，插值函数就是二元函数，即曲面。

6.1插值节点为网格节点

（1）z=interp2(x0,y0,z0,x,y,'method')

（2）pp=csape({x0,y0},z0,conds,valconds); %三次样条插值

z=fnval(pp,{x,y})

6.2 插值节点为散乱节点