MATLAB 牛顿插值法例题与程序

题目一：多项式插值某气象观测站在8：00（AM）开始每隔10分钟对天气作如下观测，用三次多项式插值函数（Newton）逼近如下曲线，插值节点数据如上表，并求出9点30分该地区的温度（x=10）。x12345678y22.523.324.421.7025.228.524.825.4二、数学原理假设有n+1个不同的节点及函数在节点上的值（x，y），（x，y），插值多项式有如下形式： （1）其中系数（i=0,1,2n）为特定系数，可由插值样条（i=0,1,2n）确定。根据均差的定义，把x看成a,b上的一点,可得f(x)= f（）+f（）fx, = f+fx, （）fx, ，x= fx, ，x+ fx, ，x(x-x)综合以上式子，把后一式代入前一式，可得到： f(x)= f+f（）+ f（）（）+ fx, ，x（）(x-x)+ fx, ，x,= N（x）+其中N（x）= f+f（）+ f（）（）+ fx, ，x（）(x-x) （2）= f(x)- N（x）= fx, ，x, （3） =（）(x-x)Newton插值的系数（i=0,1,2n）可以用差商表示。一般有 （k=0,1,2，n ） （4）把（4）代入（1）得到满足插值条件N（i=0,1,2，n）的n次Newton插值多项式N（x）=f（）+f（）+f（）（）+f（）（）（）.其中插值余项为： 介于之间。三、程序设计function y,A,C,L=newdscg(X,Y,x,M) % y为对应x的值，A为差商表，C为多项式系数，L为多项式% X为给定节点，Y为节点值，x为待求节点n=length(X); m=length(x); % n为X的长度for t=1:m z=x(t); A=zeros(n,n);A(:,1)=Y;s=0.0; p=1.0; q1=1.0; c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)- A(i-1,j-1)/(X(i)-X(i-j+1); end q1=abs(q1*(z-X(j-1);c1=c1*j; end C=A(n,n);q1=abs(q1*(z-X(n);for k=(n-1):-1:1C=conv(C,poly(X(k); d=length(C);C(d)=C(d)+A(k,k);end y(k)= polyval(C, z); %输出y值endL(k,:)=poly2sym(C); %输出多项式 syms M,X=1,3,5,7;Y=22.5,24.4,25.2,24.8;x=10; y,A,C,L=newdscg(X,Y,x,M)y = 21.7313A = 22.5000 0 0 0 24.4000 0.9500 0 0 25.2000 0.4000 -0.1375 0 24.8000 -0.2000 -0.1500 -0.0021C = -0.0021 -0.1187 1.4521 21.1688 L = - x3/480 - (19*x2)/160 + (697*x)/480 + 3387/1604、 结果分析和讨论 对于不超过三次的插值多项式，x如果选取1,3,5,7这三个点能够得到较好的三次插值多项式L=-0.0021x3-0.1187x2+1.4521x+21.1688。当x=10时，也即9点30分时的温度为21.7317度，结果分析知此值应是偏小的。对于选取不同的插值节点，能够得到不同的插值多项式，误差也不尽相同。5、 完成题目的体会与收获 对于牛顿插值法有了