606 分段lagrange三点插值

一：功能

给出n个节点X1,X2....Xn及相应的函数值Yi=f(Xi),i=1,2...n,

对p个插值点 Xj(j=1,2...p),选取最靠近插值点相邻的三个点，用lagrange三点插值公式对一元函数进行成组插值。节点等距与不等距均可。

二：算法简介

设给出n个节点，其中X1<X2<...<Xn,相应的函数值Y1,Y2...Yn,用lagrange三点插值公式，有：

L(x)=Σ(Π(X-Xj)/(Xk-Xj))Yk 注释：求和公式k由I到I+2

求积公式j由I到I+2,j≠k

其中如果X<(Xj+Xj+1)/2 那么I=j-1;

X>(Xn-2+Xn-1)/2 那么I=n-2

三：程序使用说明

1)输入参数

n为节点个数

p为插值节点个数

X(n),Y(n),一维数组，分别存放相应的节点和函数值

C(p),一维数组,存放插值点

2)输出参数

n,p分别为节点个数和插值个数

Xi,Yi节点与节点函数值

Xj,Yj,j=1,2...p,插值点与插值

BASIC源程序如下

40 input"n,p=";N,P

50 print;TAB(3);"example:";"n=";N,"p=";p

60 dim X(N),Y(N),C(P),D(P)

70 PRINT TAB(3),"JIEDIAN YU JIEDIANZHI:"

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 P

170 READ C(J)

180 NEXT J

190 PRINT

200 GOSUB 400

210 PRINT TAB(3);"CHAZHIDIAN YU CHAZHI"

220 FOR J=1 TO P

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 1.5,2.5,3.5,4.5,5.5

280 END

400 'SUBPROGRAM

410 FOR T=1 TO P

420 I=1

430 IF C(T)<.5*(X(I+1)+X(I+2)) THEN 480

440 IF C(T)>=.5*(X(N-2)+X(N-1)) THEN 470

450 I=I+1

460 GOTO 430

470 I=N-2

480 M=I+2

490 D(T)=0

500 FOR J=I TO M

510 S=1

520 FOR K=I TO M

530 IF (J-K)=0 THEN 550

540 S=S*(C(T)-X(K))/(X(J)-X(K))

550 NEXT K

560 D(T)=D(T)+S*Y(J)

570 NEXT J

580 NEXT T

590 RETURN_

#include <stdio.h> #include <stdlib.h>

void newton(float*X, float* Y, float* C, float *D, int N, int M); void main(int argc, char*argv[]) { int N , M , i; float *X, *Y, *C, *D; if(argc != 3) { do{ printf("Please Input 2 Numbers:(f.e:10 15)"); fflush(stdin); } while(scanf("%d %d", &N, &M) != 2);

} else { N = atoi(argv[1]); M = atoi(argv[2]); } printf("N = %d, M = %d\n", N, M); if((X = (float*)malloc(sizeof(float) * 2 *(N + M))) == NULL) { printf("malloc Failure!\nProgram terminate!\n"); system("pause"); exit(0); } Y = X + N; C = Y + N; D = C + M; i = 0; while(1) { if(i == N) break; printf("please input X[%d]:", i+1); if(scanf("%f" , X+i) == 1) i++; } i = 0; while(1) { if(i == N) break; printf("please input Y[%d]:", i+1); if(scanf("%f" , Y+i) == 1) i++; } i = 0; while(1) { if(i == M) break; printf("please input C[%d]:", i+1); if(scanf("%f" , C+i) == 1) i++; } for(i = 0; i < M; i++) *(D + i) = 0; newton(X, Y, C, D, N, M); for(i = 0; i < M; i++) printf("x = %0.10e y = %0.10e\n", C[i], D[i]); free(X); system("pause"); } void newton(float*X, float* Y, float* C, float *D, int N, int M) { int i, j, L; float S = 1; for(L = 1; L < N; L++) for(i = N - 1; i >= L; i--) Y[i] = (Y[i-1] - Y[i])/(X[i-L] - X[i]); for(j = 0; j < M; j++) { D[j] = Y[0]; for(i = 1; i < N; i++) { S = 1; for(L = 0; L <= i-1; L++) S *= (C[j] - X[L]); D[j] += S*Y[i]; } } }