[理学]CP030-计算物理数据插值0.ppt

1、实验数据的插值,根据这些数据，希望合理地估计出其它温度（如25摄氏度，40摄氏度）时的电阻,举例,这就是本章要讨论的“插值问题”,当精确函数 y = f(x) 非常复杂或未知时，在区间a,b上一系列节点 x0 xm处测得函数值 y0 = f(x0), , ym = f(xm)，由此构造一个简单易算的 近似函数 g(x) f(x)，满足条件 g(xj) = f(xj) (j = 0, m) (*) 这个问题称为“插值问题”,插值问题的定义,这里的 g(x) 称为f(x) 的插值函数。,节点 x0 xm称为插值节点,条件（*)称为插值条件，区间a,b称为插值区间,f(x),g(x),最常用的插值函

2、数是 ?,代数多项式,用代数多项式作插值函数的插值称为代数插值,本章主要讨论的内容,插值函数的类型,插值问题,插值法,插值函数,代数插值中的三个问题,一、插值问题解的存在唯一性？ 二、插值多项式的常用构造方法？ 三、插值函数的误差如何估计？,代数插值,代数插值问题解的存在惟一性,给定区间a,b上互异的n+1个点xjnj=0的一 组函数值f(xj)，j =0，, n，求一个n次多项式pn(x)Pn，使得 pn(xj)=f(xj)，j=0,1,，n. . (1) 令 pn(x)=a0+a1x+anxn, . (2) 只要证明Pn(x)的系数a0 ,a1, an存在唯一即可,(3),为此由插值条件(

3、1)知Pn(x)的系数满足下列n+1个代数方程构成的线性方程组 a0+a1x0+anx0n=f(x0) a0+a1x1+anx1n= f(x1) . a0+a1xn+anxnn= f(xn),代数插值问题解的存在惟一性,而ai(i=0,1,2,n)的系数行列式是Vandermonde行列式,由于xi互异，所以(4)右端不为零，从而方程组(3)的解 a0 ,a1 ,an 存在且唯一。,代数插值问题解的存在惟一性,通过解上述方程组(3)求得插值多项式pn(x)的方法并不可取.这是因为当n较大时解方程组的计算量较大，而且方程组系数矩阵的条件数一般较大（可能是病态方程组）,当阶数n越高时，病态越重。,

4、为此我们必须从其它途径来求Pn(x)：不通过求解方程组而获得插值多项式,代数插值问题解的存在惟一性,基本思想：在n次多项式空间Pn中找一组合适的基函数 0(x),1(x), n（x),使,pn(x)=a0 0(x) +a1 1(x) +an n（x),不同的基函数的选取导致不同的插值方法,Lagrange插值,Newton插值,代数插值问题解的存在惟一性,两节点-一次（线性）插值,两节点-一次（线性）插值,两节点-一次（线性）插值,两节点-一次（线性）插值,两节点-一次（线性）插值,Function y=chazhi1(x0,y0,x1,y1,x)%拉格朗日插值法 y=(x-x1)/(x0-x

5、1)*y0+( x-x0)/(x1-x0)*y1 Function y=chazhi1(x0,y0,x1,y1,x)%牛顿插值法 y=y0+(y1-y0)/(x1-x0)*( x-x0),Matlab 编程实现,三节点-二次插值,三节点-二次插值,三节点-二次插值,Function y=chazhi1(x0,y0,x1,y1,x2,y2,x)%拉格朗日插值法 y=(x-x1)*(x-x2)/(x0-x1)*(x0-x2)*y0+ (x-x0)*(x-x2)/(x1-x0)*(x1-x2)*y1+ (x-x0)*(x-x1)/(x2-x0)*(x2-x1)*y2 Function y=chazh

6、i1(x0,y0,x1,y1, x2,y2,x)%牛顿插值法 y=y0+(y1-y0)/(x1-x0)*( x-x0)+ (x-x0)*(x-x1)*(y2-y1)/(x2-x1)-(y1-y0)/(x1-x0)/(x2-x0),Matlab 编程实现,三节点-二次插值,三节点-逐次线性插值,function chazhi1 %一次插值 x0=pi/6;y0=1/2; x1=pi/4;y1=1/sqrt(2); x=pi/180*50; y=chazhi12(x0,y0,x1,y1,x) error=y-sin(x) function y=chazhi11(x0,y0,x1,y1,x)%拉格朗

7、日插值法 y=(x-x1)/(x0-x1)*y0+( x-x0)/(x1-x0)*y1; function y=chazhi12(x0,y0,x1,y1,x)%牛顿法 y=y0+(y1-y0)/(x1-x0)*( x-x0);,解:(一次插值),function chazhi2 %二次插值 x0=pi/6;y0=1/2; x1=pi/4;y1=1/sqrt(2); x2=pi/3;y2=sqrt(3)/2; x=pi/180*50; y=chazhi22(x0,y0,x1,y1,x2,y2,x)%选择插值函数 error=y-sin(x) function y=chazhi21(x0,y0,x

8、1,y1,x2,y2,x)%拉格朗日插值法 y=(x-x1)*(x-x2)/(x0-x1)*(x0-x2)*y0+(x-x0)*(x-x2)/. (x1-x0)*(x1-x2)*y1+(x-x0)*(x-x1)/(x2-x0)*(x2-x1)*y2; function y=chazhi22(x0,y0,x1,y1,x2,y2,x)%牛顿插值法 y=y0+(y1-y0)/(x1-x0)*( x-x0)+(x-x0)*(x-x1)*. (y2-y1)/(x2-x1)-(y1-y0)/(x1-x0)/(x2-x0);,解:(二次插值),N+1节点-N次插值,Lagrange插值基函数,（2）,与 节

9、点有关，而与f 无关,这里每个lj(x)都是n次多项式，且由(1)式容易验证lj(x）满足,j=0,1,，n （1）,对任意的pn(x)Pn，都有pn(x)=c0 l0(x)+c1 l1(x)+cnln(x)其中c0,c1,cn为组合系数。,可以证明函数组l0(x)，l1(x)，, ln(x) 在插值区间a,b上线性无关，所以这n+1个函数可作为Pn的一组基函数，称为Lagrange插值基函数。,Lagrange插值基函数,由Lagrange插值基函数满足(2)式可知，方程组变成,因此得到插值多项式 pn(x)= f(x0)l0(x)+f(x1) l1(x)+ f(xn) ln(x),记为Ln

10、(x) f(xj)lj(x),称Ln(x)为n次Lagrange插值多项式,插值余项 /* Remainder */,定理 若,在a , b内存在, 则在a , b上,的n+1个互异的点，对 f（x)所作的n次Lagrange插值多项式Ln (x) 有误差估计,Rolles Theorem的推论: 若 充分光滑，且,证明：由于Rn(xi) 0 ，i=0,1,，n,任意固定 x xi (i = 0, , n), 考察,(t)有 n+2 个不同的根 x0 xn x,解：,n = 1,分别利用x0, x1 以及 x1, x2 计算,利用,sin 50 = 0.7660444,利用x0, x1 作为插

11、值节点的实际误差 0.01001,利用,计算得：sin 50 0.76008,利用x1, x2作为插值节点的实际误差 0.00596,n = 2,sin 50 = 0.7660444,2次插值的实际误差 0.00061,function main %N次插值 x0=pi/6;y0=1/2; x1=pi/4;y1=1/sqrt(2); x2=pi/3;y2=sqrt(3)/2; xy=x0 x1 x2;y0 y1 y2; x=pi/180*50; y=chazhiN1(xy,x) error=y-sin(x),Matlab 编程实现,N+1节点-N次插值,function y=chazhiN1(

12、xy,x) %拉格朗日插值法 y=0; N=size(xy,2); for i=1:N L=1; for j=1:N if i=j L=L*(x-xy(1,j)/(xy(1,i)-xy(1,j); end end y=L*xy(2,i)+y; end,牛顿插值(Newtons Interpolation ),Lagrange 插值虽然易算，但若要增加一个节点时，全部基函数 li(x) 都需要重新计算。,能否重新在Pn中寻找新的基函数 ？,希望每加一个节点时，只附加一项上去即可。,1,x - x0, (x - x0)(x - x1) ，(x-x0)(x-x1) (x-xn-1)是否构成Pn的一组

13、基函数？,利用插值条件Nn(xj)=f(xj), j=0,1，n代入上式，得关于Ak (k=0,1，n)的线性代数方程组,牛顿插值法的基函数,当xj 互异时，系数矩阵非奇异，且容易求解,基函数,How complex the expression are!, 差商(亦称均差) /* divided difference */,称为在xi,xj处的1阶差商,称为在xi,xj,xk处的2阶差商,k阶差商：,利用插值条件和差商,可求出Nn(x)的系数 Ai :,Newton插值多项式,因此，每增加一个结点，Newton插值多项式只增加一项，克服了Lagrange插值的缺点。,Newton插值多项式,

14、. xk f(xk) 一阶差商 二阶差商 三阶差商 n 阶差商,差商表,例1:给定f(x)=lnx的数据表 xi 2.20 2.40 2.60 2.80 3.00 f(xi) 0.78846 0.87547 0.95551 1.02962 1.09861 1.构造差商表 2.分别写出二次、四次Newton插值多项式 解:差商表,N2(x)=0.78846 +0.43505(x-2.20) - 0.087375(x-2.20) (x-2.40） N4(x)= 0.78846 +0.43505(x-2.20） - 0.087375(x-2.20)(x-2.40) +0.0225(x-2.20)(x-2.40)(x-2.60) -0.00755(x-2.20)(x-2.40)(x-2.60)(x-2.80),function main %N次插值 x0=pi/6;y0=1/2; x1=pi/4;y1=1/sqrt(2); x2=pi/3;y2=sqrt(3)/2; xy=x0 x1 x2;y0 y1 y2; x=pi/180*50; y=chazhiN1(xy,x) error=y-sin(x),Matlab 编程实现,N+1节点-N次插值,function y=chazhiN2(xy,x)%牛顿法 N