Matlab优化设计及其应用
上QQ阅读APP看书,第一时间看更新

3.2 拉格朗日插值多项式

拉格朗日插值多项式是一种显式公式,它将pn(x)表示为一组插值基函数的线性组合。

3.2.1 线性插值

设函数y=f(x)在给定的互异节点x0,x1上的函数值分别为y0=f(x0),y1=f(x1),若能构造一个函数

p1(x)=a+bx  (3-3)

使它满足p1(x0)=y0,p1(x1)=y1,则式(3-3)所示的插值问题称为线性插值。

线性插值的几何意义是过曲线y=f(x)上的两点(x0,y0),(x1,y1)作一直线p1(x),用p1(x)来近似地替代f(x)。

对于给定的两点(x0,y0)和(x1,y1),则某一点x(x0<x<x1)处的函数值f(x)可近似表示为

l0(x)和l1(x)称为线性插值基函数。

线性插值的解可以表示为插值基函数l0(x)和l1(x)的线性组合,其组合系数为y0和y1,即

p1(x)=y0l0(x)+y1l1(x)

3.2.2 二次函数插值

二次函数插值法又称抛物线法。它的基本思路是:利用目标函数在若干点的信息和函数值,构造一个与目标函数相接近的低次插值多项式,然后求该多项式的最优解作为原函数的近似最优解。随着区间的逐次缩小,多项式的最优点与原函数最优点之间的距离逐渐缩小,直到满足一定精度要求时终止迭代。

1.二次函数插值法的基本原理

设目标函数f(x)在三点x1<x2<x3上的函数值分别为f1=f(x1),f2=f(x2),f3=f(x3),且满足f1>f2<f3,即满足函数值呈“大—小—大”的趋势。于是可通过原函数曲线上的3个点p1(x1,f1)、p2(x2,f2)、p3(x3,f3)作一条二次曲线(抛物线),如图3-4所示。此二次插值多项式为p(x)=a+bx+cx2

图3-4 二次插值函数逼近极值点的情况

为了求解p(x)的极小值,对p(x)求导数,并令其为零,即

p′(x)=b+2cx=0

解得二次函数极小点

b和c为待定参数。

根据插值条件,插值函数p(x)与原函数f(x)在插值点p1、p2、p3处函数值相等,得

求得a,b,c后即得二次插值函数极小点xp的计算公式

为便于计算,将式(3-6)改写为

式中

此时求得的xp作为函数极小点x*的近似解往往达不到精度要求,为此需缩短区间,进行多次插值计算,使xp不断逼近原函数的极小点x*。

2.二次函数插值法的迭代过程与程序框图

二次函数插值法的迭代过程如下。

(1)给定初始搜索区间[a,b]和精度ε。

(2)在区间[a,b]内取3点:x1=a,x2=0.5(a+b),x3=b,并计算函数值f1=f(x1),f2=f(x2),f3=f(x3),构成3个插值点p1(x1,f1)、p2(x2,f2)、p3(x3,f3)。

(3)按式(3-7)计算二次插值函数的极小点xp,并计算fp=f(xp)。若本步骤为第一次插值或x2点仍为初始给定点时,说明x2和xp不代表前后两次插值函数的极小点,不能进行终止判断,故进行第4步,否则转第5步。

(4)缩短搜索区间。根据原区间中x2和xp的相对位置和函数值f2和fp的比较,区间缩短有4种情况,如图3-5所示。图中阴影线部分表示舍弃的区间。

①xp>x2,f2>fp,以[x2,x3]为新区间,令x1=x2,x2=xp,x3不变,如图3-5(a);

②xp>x2,f2≤fp,以[x1,xp]为新区间,令x3=xp,x1、x2不变,如图3-5(b);

③xp≤x2,f2>fp,以[x1,x2]为新区间,令x3=x2,x2=xp,x1不变,如图3-5(c);

④xp≤x2,f2≤fp,以[xp,x3]为新区间,令x1=xp,x2、x3不变,如图3-5(d)。

对新区间的3个新点作如上的一般处理后,计算函数值f1=f(x1),f2=f(x2),f3=f(x3),返回第3步。

(5)判断是否满足精度要求。

当满足|xp-x2|≤ε时,停止迭代,把x2与xp中原函数值较小的点作为极小点;否则返回第4步,再次缩短搜索区间,直到满足精度要求为止。

图3-5 二次函数插值法区间缩短的四种情况

按上述步骤设计的二次插值函数法的程序框图如图3-6所示。根据程序框图编写的Matlab程序如下。

图3-6 二次插值法程序框图

Matlab程序如下:

function[opt_step,fo,xo]=opt_step_quad2(f,xk0,dir0,th,TolX,TolFun,MaxIter)

%opt_step_quad.m

[t012,fo,xx]=opt_range_serach1(f,xk0,dir0,th);

%search for the optimum step corresponding to minimum f(x)by quadratic approximation method

k=MaxIter;

while k>0

k=k-1;

t0=t012(1);t1=t012(2);t2=t012(3);

x0=xk0+t012(1)*dir0;x1=xk0+t012(2)*dir0;x2=xk0+t012(3)*dir0;

f0=feval(f,x0);f1=feval(f,x1);f2=feval(f,x2);

nd=[f0-f2f1-f0f2-f1]*[t1*t1t2*t2t0*t0;t1t2t0]';

t3=nd(1)/2/nd(2);x3=xk0+t3*dir0;f3=feval(f,x3);

if k<=0|abs(t3-t1)<TolX|abs(f3-f1)<TolFun

opt_step=t3;

xo=xk0+opt_step*dir0;

fo=f3;

return

else

if t3<t1

if f3<f1,t012=[t0 t3 t1];f012=[f0 f3 f1];

else t012=[t3 t1 t2];f012=[f3 f1 f2];

end

else

iff3<=f1,t012=[t1 t3 t2];f012=[f1 f3 f2];

elset012=[t0 t1 t3];f012=[f0 f1 f3];

end

end

end

end

opt_step=t3;

xo=xk0+opt_step*dir0;

fo=f3;

end

【例3-2】 用二次插值法计算目标函数

f(x)=2+x2

在区间[-2,2]内的极小点及最佳步长。

解:用户程序如下:

%opt_step_quad_test1

clc;

f=inline('2+x^2','x');

xk0=2;th=0.5;dir0=1;TolX=1e-4;TolFun=1e-4;MaxIter=50;

[opt_step,fo,xx]=opt_step_quad2(f,xk0,dir0,th,TolX,TolFun,MaxIter)

计算结果为:

opt_step=[-2]

xo=[0]

fo=[2]

3.2.3 n次拉格朗日插值多项式

1.插值多项式的存在唯一性

所谓插值多项式就是要构造一个不超过n次的多项式

pn(x)=a0+a1x+…+anxn

使其满足插值条件(3-8)。

式(3-8)是一个关于待定参数a0,a1,…,an的n+1阶线性方程组,其系数矩阵行列式为:

称为Vandermonde(范德蒙)行列式。当xi≠xj(i≠j)时,则有Vn≠0。于是,当x0,x1,…,xn为互不相同的节点时,满足条件(3-8)的插值多项式存在且唯一。

2.n次拉格朗日插值多项式

直接通过解线性方程组(3-8)来确定插值多项式的系数,一般计算工作量较大,且得出的多项式不便于应用。因此常用其他一些方法来构造插值多项式。下面我们采用所谓插值基函数的方法来构造一种具有一定特点的插值多项式。

对于给定的n+1个互异点xi(i=0,1,…,n)处的数据f(xi)=yi(i=0,1,…,n),

称li(x)为关于节点xi的n次插值基函数,且li(x)满足

再令

Ln(x)称为n次Lagrange插值多项式。可见,一个n次Lagrange插值多项式要由n+1个n次Lagrange插值基函数组合而成,其组合系数恰好是节点上的函数值。

【例3-3】 已知x0=3,x1=-1,x2=4,f(x0)=5,f(x1)=-3,f(x2)=2,求L2(x)。

解:节点基函数为

由n次Lagrange插值多项式(3-10)得

【例3-4】 根据Lagrange插值公式,用Matlab编程计算ln(x)在x=1.5的的值,差值区间x∈[1,20]。

%int_lagrange

function int_lagrange

clc;

x=1:0.5:20;

y=log(x);

xx=1.5;

n=length(x);

ff=0;

for i=1:n

p=1;

forj=1:n

if j~=i

p=p*(xx-x(j))/(x(i)-x(j));

end

end

ff=ff+p*y(i);

end

[xx,ff]

log(1.5)

计算结果为:

x=1.5000,f(x)=0.5118。