求助，有谢
603 Newton插值功能给出n个节点及相应的函数值x1,x2,…,xn;y1,y2,…,yn,
对m个插值节点xj,j=1,2,…,m,用Newton插值多项式进行成组
插值.
一. 算法简介
对已知n个节点及相应处函数值yi=f(xi),I=1,2,…,n,构
造不超过n-1次的插值多项式Φn-I(x),近似代替f(x).
当n=3时,就是构造不超过二次的插值多项式
Φ2(x)=A1+A2(x-x1)+A3(x-x1)(x-x2) (6-5)
由插值条件T2(x1)=f(x1)=y,得A1=y1.再利用插值条件
T2(x2)=y2,可得:
A2=[f(x2)-f(x1)]/(x2-x1)=f(x2,x1) (6-6)
式(6-6)右端f(x2,x1)是个记号,称为函数f(x)在x1,x2的一阶
差商,或一阶均差,再利用T2(x3)=y3,可得
A3=[f(x2)-f(x1)]/[(x3-x1)(x3-x2)]
-[(x3-x1)f(x2,x1)]/[(x3-x1)(x3-x2)]
=[f(x3,x1)-f(x1,x2)]/(x3-x1)=f(x1,x2,x3) (6-7)
式(6-7)右端f(x1,x2,x3)是个记号,称为函数f(x)在x1,x2,x3
的二阶差商或二阶均差,它是一阶均差的均差,将A1,A2,A3代入
式(6-5)即得Newton二次插值多项式.
当n>3时,我们进行推广,假定已有n-2阶差商,可递归地定
义n-1阶差商f(x1,x2,…,xn)=
[f(x1,x2,…,xn-1)-f(x1,x2,…,xn)]/(x1-xn)
由此可写出n-1阶Newton插值多项式
Φn-1=f(x1)+(x-x1)f(x1,x2)+(x-x1)(x-x2)f(x1,x2,x3)+…
(x-x1)…(x-xn-1)f(x1,x2,…,xn) (6-8)
用Φn-1(x)近似代替f(x),其误差如下:
设f(x)∈C^n-1[a,b],f^n(x)在[a,b]上存在,其中[a,b]是
包含x1,x2,…,xn在内的任意区间,则对任意给定的x∈[a,b]总
存在一点C∈(a,b)且与x有关,使
Rn-1(x)=f(x)-Tn-1(x)=f^n(C)*W(x)/n! (6-9)
其中W(x)=(x-x1)(x-x2)…(x-xn).该误差公式应用起来有一定
困难,因为实际上f(x)往往不知道.可用时候估计去估计误差.
二. 程序使用说明
1. 输入参数
n 节点个数
m 插值点个数
X(N),Y(N),一维实数组,分别存放节点及节点处的函数值
C(M) 一维实数组,存放插值点的值
数据排列顺序:x(n),y(n),c(m),分别存放在主程序250-270中
2. 输出参数
n,m分别为节点个数,插值点个数
X(N),Y(N),一维实数组,分别存放节点及节点处的函数值
C(M),D(M),一维实数组,分别存放插值点及值
当运行程序后，输入:
6 7
1.0000 2.0000 3.0000 4.0000 5.0000 6.0000
8.0000 27.0000 64.0000 125.0000 216.0000 343.0000
0 1.5 2.5 3.5 4.5 5.5 7
则运行结果为：
x=0.0000000000e+00 y=1.0000000000e+00
x=1.5000000000e+00 y=1.5625000000e+01
x=2.5000000000e+00 y=4.2875000000e+01
x=3.5000000000e+00 y=9.1125000000e+01
x=4.5000000000e+00 y=1.6637500000e+02
x=5.5000000000e+00 y=2.7462500000e+02
x=7.0000000000e+00 y=5.1200000000e+02
basic源程序:
10 '**************
20 '603 Newton *
30 '**************
40 INPUT "n,m=";N,M
50 PRINT;TAB(3);"n=";N,"m=";M
60 DIM X(60),Y(60),C(M),D(M)
70 PRINT TAB(3);
80 FOR I=1 TO N
90 READ X(I):PRINT USING" ###.####";X(I);
100 NEXT I
110 PRINT
120 FOR I=1 TO N
130 READ Y(I):PRINT USING" ###.####";Y(I);
140 NEXT I
150 PRINT
160 FOR J=1 TO M
170 READ C(J)
180 NEXT J
190 PRINT
200 GOSUB 400
210 PRINT TAB(3);
220 FOR J=1 TO M
230 PRINT "X=";C(J),"Y=";D(J)
240 NEXT J
250 DATA 1,2,3,4,5,6
260 DATA 8,27,64,125,216,343
270 DATA 0,1.5,2.5,3.5,4.5,5.5,7
280 END
400 '子程序
410 FOR L=2 TO N
420 FOR I=N TO L STEP-1
430 Y(I)=(Y(I-1)-Y(I))/(X(I-L+1)-X(I))
440 NEXT I
450 NEXT L
460 FOR J=1 TO M
470 D(J)=0
480 FOR I=1 TO N
490 S=1
500 IF I=1 THEN 540
510 FOR L=1 TO I-1
520 S=S*(C(J)-X(L))
530 NEXT L
540 D(J)=D(J)+S*Y(I)
550 NEXT I
560 NEXT J
570 RETURN
