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

实验数据的插值 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) (*) 这个问题称为“插值问题” 插值问题的定义 这里的 g(x) 称为f(x) 的插值函数。 节点 x0 xm称为插值节点, 条件（*)称为插值条件，区间a,b称为插值区间 x0x1 x2 x3x4 x f(x) g(x) 最常用的插值函数是 ? 代数多项式 用代数多项式作插值函数的插值称为代数插值 本章主要讨论的内 容 插值函数的类型 插值问题 插值法 插值函数 代数插值中的三个问题 n一、插值问题解的存在唯一性？ n二、插值多项式的常用构造方法？ n三、插值函数的误差如何估计？ 代数插值 代数插值问题解的存在惟一性 给定区间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)知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越高时，病态越重。 为此我们必须从其它途径来求 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-x1)*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=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 编程实现 三节点-二次插值 三节点-逐次线性插值 例：已知 分别利用 sin x 的1次、2次 Lagrange 插值和牛顿插 值计算 sin 50。 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)%拉格朗日插值法 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,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,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） 与 节点有关，而与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(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 例：已知 分别利用 sin x 的1次、2次 Lagrange 插值计算 sin 50， 并估计误差。 解： n = 1分别利用x0, x1 以及 x1, x2 计算 利用 sin 50 = 0.7660444 利用x0, x1 作为插值节点的实际误差 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(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的一组基函数？ 利用插值条件Nn(xj)=f(xj), j=0,1，n代入上式 ，得关于Ak (k=0,1，n)的线性代数方程组 牛顿插值法的基函数 当xj 互异时，系数矩阵非奇异，且容易求解 基函数 How complex the expression are! It is not a difficult thing for a mathematician. We can use notation 差商(亦称均差) /* divided difference */ 称为在xi,xj处的1阶差商 称为在xi,xj,xk处的2阶差商 k阶差商： 利用插值条件和差商,可求出Nn(x)的系数 Ai : Newton插值多项式 因此，每增加一个结点，Newton插值多项式只增加 一项，克服了Lagrange插值的缺点。 Newton插值多项式 . 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=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,