Lagrange插值法1.ppt

1、在生产和实验中，函数f (x)无表达式, 只知道f(x)在一些给定点的函数值(或其导数值) ，或者其表达式复杂不便于计算，此时我们希望建立一个简单而又便于计算的函数(x)，使其近似的代替f(x).,第四章 插值与拟合,插值法：插值法是利用函数f(x)在一组节点x0 , x1, . , xn上的值y0 , y1 , , yn，构造一个插值函数(x)来逼近已知函数f(x) ，并要求插值函数与已知函数在节点处的函数值相同,即:,(x0)=y0 , (x1)=y1 , , (xn)=yn,求近似函数(x)的方法一般分为两类: 一类是插值,另一类是拟合.,引例及问题综述,插值法：插值法是利用函数在一组节

2、点上的值，构造一个插值函数(x)来逼近已知函数f(x) ，并要求插值函数与已知函数在节点处的函数值相同。,曲线拟合方法:不要求近似函数(x)所表示的函数严格通过已知数据点(xi,yi),而是通过观察这些点的分布规律，选择某种能描述这一近似规律的函数(x)来逼近函数f(x)，然后按照某种原则使逼近效果总体上尽可能好，其中最常见的原则就是最小二乘原理。,但在实际应用中，节点上的函数值通常不是精确值，而是由实验或测量得到的数据，不可避免的带有误差，如果用插值法，会保留这些误差，影响精度。为了尽量减少这种测量误差，人们又提出了另外一种构造近似函数的方法曲线拟合方法。,原始数据的散点图,分段线性插值,二

3、次分段拟合曲线,已知数据点(xi,yi),艾滋病疗法的评价及疗效的预测 2006高教社杯全国大学生数学建模竞赛题目,艾滋病的医学全名为“获得性免疫缺损综合症”，英文简称AIDS，它是由艾滋病毒（医学全名为“人体免疫缺损病毒”, 英文简称HIV）引起的。这种病毒破坏人的免疫系统，使人体丧失抵抗各种疾病的能力，从而严重危害人的生命。人类免疫系统的CD4细胞在抵御HIV的入侵中起着重要作用，当CD4被HIV感染而裂解时，其数量会急剧减少，HIV将迅速增加，导致AIDS发作。 艾滋病治疗的目的，是尽量减少人体内HIV的数量，同时产生更多的CD4，至少要有效地降低CD4减少的速度，以提高人体免疫能力。

4、迄今为止人类还没有找到能根治AIDS的疗法，目前的一些AIDS疗法不仅对人体有副作用，而且成本也很高。许多国家和医疗组织都在积极试验、寻找更好的AIDS疗法。,附件1: 同时服用3种药物（zidovudine(齐多夫定 ), lamivudine(拉米夫定 ),indinavir(茚地那韦 )）的300多名病人每隔几周测试的CD4和HIV的浓度。 第1列是病人编号，第2列是测试CD4的时刻（周），第3列是测得的CD4（乘以0.2个/ml），第4列是测试HIV的时刻（周），第5列是测得的HIV（单位不详）。,现在得到了美国艾滋病医疗试验机构ACTG公布的两组数据。(略),PtID CD4Date

5、CD4CountRNADateVLoad 234240 178 0 5.5 234244 228 4 3.9 234248 126 8 4.7 2342425 171 25 4 2342440 99 40 5 0 14 0 5.3 ,附件2: 1300多名病人按照4种疗法服药大约每隔8周测试的CD4浓度。 第1列是病人编号，第2列是4种疗法的代码： 1 = 600mg zidovudine(齐多夫定) 与400mg didanosine(去羟肌苷)按月轮换使用； 2 = 600mg zidovudine 加2.25mg zalcitabine(双脱氧胞苷 )； 3 = 600mg zidovu

6、dine 加400mg didanosine； 4 = 600mg zidovudine 加400mg didanosine 加400mg nevirapine(奈韦拉平 )。 第3列是病人年龄，第4列是测试CD4的时刻（周），第5列是测得的CD4，取值log(CD4+1).,ID 疗法 年龄 时间 Log(CD4 count+1) 1236.42710 3.1355 1236.42717.57143.0445 1236.427115.57142.7726 1236.427123.57142.8332 1236.427132.57143.2189 .,请你完成以下问题： （1）利用附件1的数据

7、，预测继续治疗的效果，或者确定最佳治疗终止时间（继续治疗指在测试终止后继续服药，如果认为继续服药效果不好，则可选择提前终止治疗）。 （2）利用附件2的数据，评价4种疗法的优劣（仅以CD4为标准），并对较优的疗法预测继续治疗的效果，或者确定最佳治疗终止时间。 （3） 艾滋病药品的主要供给商对不发达国家提供的药品价格如下：600mg zidovudine 1.60美元，400mg didanosine 0.85美元，2.25 mg zalcitabine 1.85美元，400 mg nevirapine 1.20美元。如果病人需要考虑4种疗法的费用，对（2）中的评价和预测（或者提前终止）有什么改变

8、。,生长阻滞模型,缺少数据,缺少数据是用样条插值函数求出来的,高等教育学费问题探讨,*数据散点图 O 拟合曲线图,(*其中部分数据是用样条插值函数求出来的),样条插值,插值法是函数逼近的重要方法之一，有着广泛的应用 。 基本思想:就是利用函数f (x)在一些给定点的函数值(或其导数值) ，建立一个简单而又便于计算的函数(x)，使其近似的代替f(x) . 插值法有很多种,其中以拉格朗日(Lagrange)插值和牛顿(Newton)插值为代表的多项式插值最有特点,常用的插值还有Hermit插值,分段插值和样条插值.,插值法,插值法 设函数 y=f(x) 在区间a,b上有定义,且已知f(x)在a,b

9、上n+1个互异点x0 , x1, . , xn 处的值 y0 , y1 , , yn , 若存在一个近似函数 (x) ，满足 ： (x0)=y0 , (x1)=y1 , , (xn)=yn , (*) 则称(x)为f(x)的插值函数，f(x) 称为被插值函数，x0 , x1, . , xn 称为插值节点。(*)式称为插值条件。而误差函数R(x)=f(x) - (x)称为插值余项。,基本概念,满足同一插值条件的插值函数(x)有许多类型，如：多项式函数类型、三角函数类型、指数函数类型等，常用的插值函数是多项式,我们称其为代数插值（或多项式插值）。,(x)作为函数 y=f(x) 的近似表达式(近似函

10、数),满足: y= f(x) (x)我们把构造满足插值条件的近似函数(x)，称为插值问题。,我们这章只讨论代数插值：插值多项式Pn(x)是否存在?如何求解?插值误差和余项如何估计?,最简单的插值函数是代数多项式 Pn(x)=a0+a1x+anxn, . (1) 这时插值问题变为:求一个n次多项式Pn(x),使其满足插值条件: pn(xi)=yi, i= 0,1,2,，n, (2),pn(x0)=y0,pn(x1)=y1,pn(xn)=yn,而ai(i=0,1,2,n)的系数行列式是Vandermonde行列式,由于xi互异，所以(4)右端不为零，从而方程组(3)的解 a0 ,a1 ,an 存在

11、且唯一。解出ai(i=0,1,2,n), Pn(x)就可构造出来了。但遗憾的是方程组(3)是病态方程组,当阶数n越高时，病态越重。为此我们从另一途径来寻求获得Pn(x) 的方法-Lagrange插值和Newton插值。,Lagrange插值,1 Lagrange插值多项式,设函数y=f(x)在区间a,b上有定义，且已知在点ax0x1.xnb上的函数值为y0,y1,y2,.,yn，求一个次数不超过n的多项式 Ln(x)=a0+a1x+.+anxn (1)* 使得 Ln(xi)=yi (i=0,1,2,.,n) (2)* 成立。称(1)*式为满足插值条件(2)*的拉格朗日插值多项式。,由两点式,可

12、求L1(x)的表达式,整理得:,另一种推导方式:令L1(x)=l0(x)y0+l1(x)y1,一、拉格朗日插值多项式的构造,线性插值多项式,1、线性插值,先从最简单的线性插值(n=1)开始.即有：a=x0，b=x1, 这时插值问题就是求一次多项式L1(x)(=c+dx )，使其满足条件： L1(x0)=y0 , L1(x1)=y1 ,即l0(x)含有因子x-x1， l1(x)含有因子(x-x0)，,l0(x), l1(x)称为以x0 , x1 为节点的线性插值基函数,这样得到一次插值多项式:,令L1(x)=l0(x)y0+l1(x)y1,l0(x0)=1, l1(x0)=0, l0(x1)=0

13、, l1(x1)=1.,令 l0(x)=0 (x-x1) , l1(x)=1 (x-x0)，利用 l0(x0)=1 和 l1(x1)=1确定其中的系数0， 1得：,线性插值多项式,0=1/(x0-x1),故有:,由L1(x0)=y0 , L1(x1)=y1得,1=1/(x1-x0),令 L2(x)=l0(x)y0 + l1(x)y1 + l2(x)y2 。,线性插值仅仅用两个节点上的信息，精确度较差。为了提高精确度，我们进一步考察以下三点的插值问题：,求二次多项式 L2(x) =a0 + a1x + a2x2，使其满足条件：L2(x0)=y0 , L2(x1)=y1 , L2(x2)=y2,类

14、似的可以得出 l1(x) , l2(x) :,这样 l0(x)含有 x-x1 , x-x2 两个因子，令 l0(x)=0 (x-x1)(x-x2) ,利用 l0(x0)=1 确定其中的系数0，得,2、抛物线插值,l0 (x0 )=1 , l1 (x0 )=0 , l2 (x0 )=0 , l0 (x1 )=0 , l1 (x1 )=1 , l2 (x1 )=0 , l0 (x2 )=0 , l1 (x2)=0 , l2 (x2 )=1 .,于是 (x-x1)(x-x2) (x-x0)(x-x2) (x-x0)(x-x1) L2(x)=-y0 + -y1 + -y2 .(4) (x0-x1)(x

15、0-x2) (x1-x0)(x1-x2) (x2-x0)(x2-x1) l0(x) , l1(x) , l2(x) 称 为以 x0 , x1 , x2为节点的抛物 线插值基函数。 3、n次多项式插值：仿照线性插值和二次插值的办法， 进一步讨论一般形式的 n 次多项式 Ln(x)=a0 +a1x +a2x2 + + anxn , 使其满足 Ln(x0 )=y0 , Ln(x1 )=y1 , . , Ln(xn )=yn .(5) 我们仍从构造插值基函数着手，先对某个固定的下标 j ，作 n 次多项式 l j(x) , 使其满足条件 可求得,二次插值多项式,.(*) 公式（*）就是Lagrange

16、插值多项式，lj(x)称为以x0 , x1 ,. , xn为节点的Lagrange插值基函数。,将lj(x)代入:,例1 求过点(2,0)(4,3)(6,5)(8,4)(10,1)的拉格朗日型插值多项式,解: 用5个节点作4次插值多项式,例1 求过点(2,0)(4,3)(6,5)(8,4)(10,1)的拉格朗日型插值多项式,解: 用5个节点作4次插值多项式,L4(x)= y0l0(x)+ y1l1(x)+ y2l2(x)+ y3l3(x)+y4l4(x),例2: 已给sin0.32=0.314567,sin0.34=0.333487, sin0.36=0.352274, 用线性插值及抛物插值计

17、算 sin0.3367 的值。,y1 - y0 sin0.3367L1(0.3367)=y0+(0.3367 -x0) x1 - x0 0.01892 =0.314567+ (0.0167) =0.330365 . 0.02,用线性插值计算，取 x0=0.32 及 x1=0.34 ,解: 由题意取,得:,由线性插值公式,用抛物插值计算 sin0.3367时，因为,这个结果与六位有效数字的正弦函数表完全一样，这说明查表时用二次插值精度已相当高了。插值多项式的存在性?其截断误差是多少?,所以,定理1 假设x0 ,x1,xn 是n+1个互异节点,函数f(x)在这组节点的值f(xj)(j=0,1,n)

18、是给定的，那么存在唯一的次数不超过n 的多项式Pn (x)满足 Pn (xj)=f(xj), j=0,1,n,二 、 插值多项式的存在性和唯一性,证明：由插值函数的构造知n次多项式Pn (x)的存在性，由代数基本定理可证明它的唯一性。,函数f(x)的n次插值多项式Ln(x)只是在节点处有:,Ln (xj)=f (xj), j=0,1,n,若xxj,一般Ln (x) f (x), 令 Rn(x)=f (x)-Ln(x),称Rn(x)为插值多项式的余项,三 、 Lagrange插值的插值余项 (截断误差) 定理2：设Ln(x)是过点x0 ，x1 ，x2 ，xn的 n 次插值多项式， f (x)在a

19、，b上有n阶连续导数，在(a，b)内存在n+1阶导数，其中a，b是包含点x0 , x1 , x2 , xn的任一区间，则对任意给定的xa，b，总存在一点(a，b)（依赖于x）使 其中 。,证明: Rn(x) = f(x) - Ln(x), 当x=xi时，显然有： Rn(xi) =f(xi) - Ln(xi), =0 ， n+1 (xi)=0 （ i=0,1,n）结论成立,当xxi时，Rn(xi ) =0 ，n+1 (xi)=0 （ i=0,1,n） 可设Rn(x)=K(x) n+1 (x),需证:K(x)=f(n+1)( )/(n+1)!,现在a,b上任意固定一点x,引进辅助函数 g(t)=f

20、(t)- Ln(t)-K(x) n+1 (t), (*) 则g(t)在a,b上具有n阶连续导数，在(a,b)内存在n+1阶导数，当 t= x, x0, x1, xn 时，g(t)=0,，即g(t)在(a,b)内有n+2个零点,由Rolle定理知g(t)在(a,b)内有n+1个零点,如此反复，最后可推知g(n+1)(t)在(a,b)内有1个零点, 即有g(n+1)( )=0, a b。这样，由(*)式便有,现在a,b上任意固定一点x,引进辅助函数 g(t)=f(t)- Ln(t)-K(x) n+1 (t), (*) 有g(n+1)( )=0, a b。这样，由(*)式便有,因为 n+1 (t)是

21、n+1次多项式， n+1 (n+1)(t)=(n+1)!,又因为Ln(t)是次数为n的多项式,因此Ln (n+1)(t) = 0 。 由此得： K(x)=f(n+1)( )/(n+1)! . 代入Rn(x)=K(x) n+1(x)，定理得证。,需证:K(x)=f(n+1)( )/(n+1)!,上式称为带余项的Lagrange插值公式。,当f(x)具有n+1阶导数，才有余项表达式成立，且余项为：,由此得：,那么插值多项式Ln(x)逼近f(x)的截断误差是,注：1）余项表达式只有在 f(x) 的n+1阶导数存在时才能应用。,2） 在 （a，b）内的具体位置通常不可能给出，如果我们可以求出,例2:

22、已给sin0.32=0.314567,sin0.34=0.333487, sin0.36=0.352274, 用线性插值及抛物插值计算 sin0.3367的值并估计截断误差。,解： 用线性插值计算，取 x0=0.32 及 x1=0.34 , 又由公式得 y1 - y0 sin0.3367L1(0.3367)=y0+(0.3367 -x0) x1 - x0 0.01892 =0.314567+ (0.0167) =0.330365 . 0.02,其截断误差得,其中,其截断误差得,R1(0.3367)=sin 0.3367 L1(0.3367)1/2(0.3335)(0.3367-0.32)(0.

23、 3367-0.34) 1/2(0.3335)(0.0167)(0.0033)0.92105,用线性插值计算，取 x0=0.32 及 x1=0.34 ,其中 ，,可取，,因 f(x)=sinx，f (x)= -sinx，,于是,sin0.3367,若取x1=0.34，x2=0.36为节点，则线性插值为,其截断误差为,其中,于是,注：从计算可看出，插值误差与节点xi和点x之间的距离有关，节点距离x越近，一般插值误差小。当 x位于插值区间a,b以外（称为外插）一般会引起较大误差。应避免外插。,R1(0.3367)=sin 0.3367 L1(0.3367)1/2(0.3335)(0.0167)(0

24、.0033)0.92105,误差值比较：,取 x0=0.32 及 x1=0.34,取x1=0.34，x2=0.36为节点,其截断误差为,用抛物插值计算 sin0.3367时，因为,所以,其中,其截断误差为,其中,于是,显然，抛物插值比线性插值误差小的多。,线性插值误差:R1(0.3367)0.92105,练习: 已测得某地大气压强随高度变化的一组数据,高度(m) 0 100 300 1000 1500 2000 . 压强 (kgf/m2) 0.9689 0.9322 0.8969 0.8515 0.7984 0.7485 试用二次插值法求1200米处的压强值.,作业：P107 1、2,解：设x

25、为高度，y为大气压强的值， 选取(1000,0.8515) ,(1500,0.7984), (2000,0.7485)三点构造二次插值多项式 (x-x1)(x-x2) (x-x0)(x-x2) (x-x0)(x- x1) L2(x)=- - y0 + - y1 + - y2 (x0- x1)(x0-x2) (x1 -x0)(x1 -x2) (x2-x0)(x2- x1) 代入已知的数值，得 L2(1200)=0.8515(1200-1500)(1200- 2000)/(1000-1500)(1000-2000)+0.7984(1200-1000)(1200-2000)+0.7485(1200-

26、1000)(1200-1500)/(2000-1000)(2000-1500)=300*800*0.8515/500/1000+200*800*0.7984/500/500-200*300*0.7485/500/1000=0.82980 所以 y(1200) L2(1200)= 0.82980 (kgf/m2),T=T*(x-xk)/(xj-xk) k=0,1,2,.,j-1, j+1,.,n,j=0,1,2,n,T=1,j=1,内层循环,if k=j,用for循环语句(对k),程序设计,对每一个T=T(x-xk)/(xj-xk) k=0,1,2,.,j-1, j+1,.,n,j=0,1,2,

27、n,外层循环,L=T*y0(j)+L,最后输出L,用for循环语句,L=0,y0(j)=yj,开始,输入xi ,yi i =0,1,2,.,n,L=0,j=0,T=1,T=T(x-xk)/(xj-xk) k=0,1,2,.,j-1, j+1,.,n,L=L+T*yj,j=n,j=j+1,输出L,结束,N,Y,编程思想,返回首页,function y=lagrange(x0,y0,x) n=length(x0);m=length(x); for i=1:m z=x(i); L=0.0; for j=1:n T=1.0; for k=1:n if k=j T=T*(z-x0(k)/(x0(j)-x

28、0(k); end end L=T*y0(j)+L; end y(i)=L; end,MatLab实现,% x0为全部插值节点(向量) % y0为插值节点处的函数值(向量) % x为被估计函数自变量(可为向量) % y为被估计函数x处的值(可为向量),% 输出被估计函数在x(i)处的值,x0=0.4: 0.1: 0.8; y0=-0.916291 -0.693147 -0.510826 -0.356675 -0.223144; lagrange(x0,y0,0.54),在MatLab命令窗口输入：,例 给出f (x)=lnx的数值表，用Lagrange插值计算 ln(0.54)的近似值。,输出

29、结果: -0.616143,精确解: -0.616186,( -0.6161),x0=0.4: 0.1: 0.8; y0=-0.916291 -0.693147 -0.510826 -0.356675 -0.223144; lagrange(x0,y0,0.54,0.65) ans = -0.6161 -0.4308, lagrange(x0,y0,0.54,0.65,0.78) ans = -0.6161 -0.4308 -0.2484,实验4.1（观察龙格(Runge)现象实验） 实验目的：观察拉格朗日插值的龙格(Runge)现象. 实验内容： 1、给出拉格朗日插值多项式的算法流程和相关程序； 2、对于函数 进行拉格朗日插值，取不同的节点数，在区间-5,5上取等距间隔的节点为插值点，把f(x)和插值多项式的曲线画在同一张图上进行比较。(a可以取任意值)具体步骤： 1） a=1时， i）取n=4，作出f(x)和插值多项式的曲线图； ii）取n=10，作出f(x)