CP030-计算物理数据插值课件.ppt

1、实验数据的插值n根据这些数据，希望合理地估计出其它温度（如25摄氏度，40摄氏度）时的电阻举例这就是本章要讨论的“插值问题”n已经测得在某热敏电阻在不同温度下的电阻如下：电阻（欧姆）600 300 100 50 10 5 温度（摄氏度）10 30 50 70 90 100 当精确函数当精确函数 y=f(x)非常复杂或未知时，在区非常复杂或未知时，在区间间a,b上一系列节点上一系列节点 x0 xm处测得函数值处测得函数值 y0=f(x0),ym=f(xm)，由此构造一个简单易算的，由此构造一个简单易算的近似函数近似函数 g(x)f(x)，满足条件，满足条件 g(xj)=f(xj)(j=0,m)(

2、*)这个问题称为这个问题称为“插值问题插值问题”插值问题的定义这里的这里的 g(x)称为称为f(x)的的插值函数插值函数。节点节点 x0 xm称为称为插值节点插值节点,条件（条件（*)称为称为插值条件插值条件，区间，区间a,b称为称为插值区间插值区间x0 x1x2x3x4 xf(x)g(x)最常用的插值函数是最常用的插值函数是?代数多项式代数多项式用代数多项式作插值函数的插值称为用代数多项式作插值函数的插值称为代数插值代数插值本章主要讨论的内本章主要讨论的内容容插值函数的类型插值函数的类型插值问题插值法插值函数代数插值中的三个问题n一、插值问题解的存在唯一性？n二、插值多项式的常用构造方法？n

3、三、插值函数的误差如何估计？代数插值代数插值代数插值问题解的存在惟一性代数插值问题解的存在惟一性 给定区间给定区间a,b上互异的上互异的n+1个点个点xj的一的一 组函数值组函数值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)为此由插值条件为此由插值条件(1)(1)知知P Pn n(x)(x)的系数满足下列的系数满足下列n+1n+1个个代数方程构成的线性方程组代数方程构成的线性方程组

4、a0+a1x0+anx0n=f(x0)a0+a1x1+anx1n=f(x1).a0+a1xn+anxnn=f(xn)代数插值问题解的存在惟一性代数插值问题解的存在惟一性200021110121.1.V(,.,).1.nnnnnnnxxxxxxx xxxxx110()niijijxx而而ai(i=0,1,2,n)的系数行列式是的系数行列式是Vandermonde行列式行列式由于xi互异，所以(4)右端不为零，从而方程组(3)的解 a0,a1,an 存在且唯一。代数插值问题解的存在惟一性代数插值问题解的存在惟一性 通过解上述方程组通过解上述方程组(3)求得插值多项式求得插值多项式pn(x)的方法的

5、方法并不可取并不可取.这是因为当这是因为当n较大时解方程组的计算量较大时解方程组的计算量较大，而且方程组系数矩阵的条件数一般较大较大，而且方程组系数矩阵的条件数一般较大（可能（可能是病态方程组）是病态方程组）,当阶数当阶数n越高时，病态越越高时，病态越重。重。为此我们必须从其它途径来求为此我们必须从其它途径来求Pn(x)：不通过求解方程组而获得不通过求解方程组而获得插值多项式插值多项式代数多项式插值存在的问题代数多项式插值存在的问题基本思想：基本思想：在n次多项式空间P Pn n中找一组合适的基函数中找一组合适的基函数 0 0(x),(x),1 1(x),(x),n n（x),x),使使pn(

6、x)=a0 0 0(x)(x)+a1 1 1(x)(x)+an n n（x)x)不同的基函数的选取导致不同的不同的基函数的选取导致不同的插值方法插值方法Lagrange插值插值Newton插值插值代数多项式插值的基本思想代数多项式插值的基本思想两节点两节点-一次（线性）插值一次（线性）插值两节点两节点-一次（线性）插值一次（线性）插值两节点两节点-一次（线性）插值一次（线性）插值两节点两节点-一次（线性）插值一次（线性）插值两节点两节点-一次（线性）插值一次（线性）插值Function y=chazhi1(x0,y0,x1,y1,x)%拉格朗日插值法拉格朗日插值法y=(x-x1)/(x0-x1

7、)*y0+(x-x0)/(x1-x0)*y1Function 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

8、-x1)*y2Function y=chazhi1(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 编程实现编程实现三节点三节点-二次插值二次插值三节点三节点-逐次线性插值逐次线性插值例：已知例：已知233sin,214sin,216sin分别利用分别利用 sin x 的的1次、次、2次次 Lagrange 插值和牛顿插插值和牛顿插值计算值计算 sin 50。function chazhi1%一次插值一次插

9、值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)%拉格朗日插值法拉格朗日插值法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=

10、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,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=chazhi22(x0,y0,x1,y1,x2

11、,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+1N+1节点节点-N-N次插值次插值Lagrange插值基函数插值基函数011101110()()()()()()()()()()()jjnjjjjjjjjnniijiijxxxxxxxxxxl xxxxxxxxxxxxxxxjijixlij,1,0)(（2）与与 节点节点有关，而与有关，而与f 无关无关这里每个这里每个lj(x)都是都是n n次多项式，且由次多项式，

12、且由(1)(1)式容易验证式容易验证lj(x）满足满足j=0,1,，n （1）对任意的对任意的pn(x)Pn，都有都有pn(x)=c0 l0(x)+c1 l1(x)+cnln(x)其其中中c0,c1,cn为组合系数。为组合系数。0010000011111101()()()()()()()()()()()()nnnnnnnnl xl xl xcf xl xl xl xcf xl xl xl xcf x 可以证明函数组可以证明函数组l0(x)，l1(x)，,ln(x)在插值区间在插值区间a,b上上线性无关，所以这线性无关，所以这n+1个函数可作为个函数可作为Pn的一组基函数，称为的一组基函数，称为

13、Lagrange插值基函数。插值基函数。Lagrange插值基函数插值基函数0011()100()010()001nncfxcfxcfx 由由LagrangeLagrange插值基函数满足插值基函数满足(2)(2)式可知，方程组变成式可知，方程组变成因此得到插值多项式因此得到插值多项式pn(x)=f(x0)l0(x)+f(x1)l1(x)+f(xn)ln(x)记为记为Ln(x)f(xf(xj j)l)lj j(x)(x)称称L Ln n(x)(x)为为n n次次LagrangeLagrange插值多项式插值多项式插值余项插值余项/*Remainder*/定理定理 若若(1)()nfx在在a,b

14、内存在内存在,则在则在a,b上上(1)0()()()()()(1)!nnnniifR xf xL xxxn的的n+1个互异的点，对个互异的点，对 f（x)所作的所作的n次次Lagrange插插值多项式值多项式Ln(x)有误差估计有误差估计 Rolles Theorem的推论的推论:若若 充分光滑，且充分光滑，且0)()(0 nxx )(x证明：由于证明：由于Rn(xi)0，i=0,1,=0,1,，n n()()(0nRxxu xxnii任意固定任意固定 x xi (i=0,n),考察考察00()()()()()()()()nnniniiitR tu xt xf tL tu xt x(t)有有

15、n+2 个不同的根个不同的根 x0 xn x(1)()0,(,)nxxa b(1)()()(1)!nfu xn例：已知例：已知233sin,214sin,216sin分别利用分别利用 sin x 的的1次、次、2次次 Lagrange 插值计算插值计算 sin 50，并估计误差。并估计误差。解：解：n=1分别利用分别利用x0,x1 以及以及 x1,x2 计算计算4,610 xx利用利用216/4/6/214/6/4/)(1 xxxL015sin50()0.7761418L(2)1()13()()(),sin2!6422xxfR xxx150.01319()0.0076218R sin 50 =

16、0.7660444利用利用x0,x1 作为插值节点的实际误差作为插值节点的实际误差 0.010010.010013,421 xx利用利用 计算得：计算得：sin 50 0.76008,150.005380.0066018R利用利用x1,x2作为插值节点的实际误差作为插值节点的实际误差 0.00596 0.00596n=24363264634643643634()()()()11()()()2()()2()()3()()2xxxxL xxx23cos21;)3)(4)(6(!3cos)(2 xxxxxxR 250.000440.0007718Rsin 50 =0.76604442次插值的实际误差

17、次插值的实际误差 0.00061 0.00061025sin50()0.7654318Lfunction 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+1N+1节点节点-N-N次插值次插值 function y=chazhiN1(xy,x)%拉格朗日插值法拉格朗日插值法y=0;N=size(xy,2);for i=1:N L=1;for j=1:

18、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,1,x-x0,(x-x0)(x-x1)，(x-x0)(x-x1)(x-xn-1)是否构成是否构成P Pn n的

19、一组基函数的一组基函数？01020101()()()().().()nnnN xA A x xA x x x xA x xx x 利用插值条件利用插值条件Nn(xj)=f(xj),j=0,1，n代入上式，代入上式，得关于得关于Ak(k=0,1，n)的线性代数方程组的线性代数方程组牛顿插值法的基函数牛顿插值法的基函数0010111000100()10()1()()nninniAf xxxAf xxxxxAf x当当xj 互异时，系数矩阵非奇异，且容易求解互异时，系数矩阵非奇异，且容易求解基函数基函数10110()()f xf xAxx10212202110()()()()()/()f xf xf

20、 xf xAxxxxxxHow complex the expression are!It is not a difficult thing for a mathematician.We can use notation00()Af x 差商差商(亦称均差亦称均差)/*divided difference*/),()()(,jijijijixxjixxxfxfxxf 称为在称为在xi,xj处的处的1阶差商阶差商)(,kixxxxfxxfxxxfkikjjikji 称为在称为在xi,xj,xk处的处的2阶差商阶差商k阶差商：阶差商：01112010,kkkkf xxxf x xxf xxxxx利

21、用插值条件和差商,可求出Nn(x)的系数 Ai:00100011()(),(),()()()nnnN xf xf x x x xf xx x x x xx x 00010101(),nnAf xf xAf x xAf x xxNewton插值多项式0010010001010101()()()(),()(),()()(),()()()kkkkN xf xN xf xf x x x xN xf x x x xNxN xf xxx x x xx x 因此，每增加一个结点，Newton插值多项式只增加一项，克服了Lagrange插值的缺点。Newton插值多项式.xk f(xk)一阶差商一阶差商 二阶

22、差商二阶差商 三阶差商三阶差商 n 阶差商阶差商差商表差商表0123nxxxxx0123 nf xf xf xf xf x0112231,nnf x xf x xf x xf xx01212321,nnnf x x xf x x xf xxx0123321,nnnnf x x x xf xxxx0 1,nf x xx例例1:1:给定f(x)=lnx的数据表xi 2.20 2.40 2.60 2.80 3.00f(xi)0.78846 0.87547 0.95551 1.02962 1.098611.构造差商表2.分别写出二次、四次Newton插值多项式解解:差商表差商表00755.001646

23、.006400.034495.009861.100.302250.0073875.037055.002962.180.2087375.040010.095551.060.243505.087547.040.278846.020.2四阶差商三阶差商二阶差商一阶差商iixfxN2(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)

24、(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+1N+1节点节点-N-N次插值次插值function y=chazhiN2(xy,x)%牛顿法牛顿法N=size(xy,2);f=zeros(N,N);%差商矩阵差商矩阵f(:,1)=xy(2,:);x0=xy(1,:);y=f(1);sx=1;for i=2:N sx=sx*(x-x0(i-1);for j=i:N f(j,i)=(f(j,i-1)-f(j-1,i-1)/(x0(j)-x0(j-i+1);end y=y+f(i,i)*sx;end其他插值方法n样条插值n埃尔米特（Hermitte）插值本章作业n给定某系统的运动满足函数y=exp(-t)*sint,求该函数的拉格朗日和牛顿插值函数。n要求：1、根据该函数的图形特征，设定一些节点。2、根据节点个数，写出拉格朗日和牛顿插值函数的程序。3、在同一坐标中绘出被插函数、拉格朗日插值函数、牛顿插值函数的曲线。4、分析其插值效果及影响因素。