用C语言编写一个线性插值程序

#include <stdio.h>

double Lerp(double x0,double y0,double x1,double y1,double x)

{

    double dy = y1 - y0

    if(dy == 0){

        printf("除0错误！\n")

        return 0

    }

    return x * (x1 - x0) / dy

}

int main()

{

    double x0,x1,y1,y0,x,y

    printf("Inptu x0 y0 x1 y1 x:")

    scanf("%lf %lf %lf %lf %lf",&x0,&y0,&x1,&y1,&x)

    y = Lerp(x0,y0,x1,y1,x)

    printf("y = %lf\n",y)

    return 0

}

这是一元全区间等距插值子程序，X1和H为实数，分别为等距节点中的第一个节点值和等距节点的步长；N为整数，等距节点的个数；Y(N)，存放N个等距节点上的函数值；T是指定插值点的值；Z返回指定插值点T处的函数近似值。

SUBROUTINE EELGR(X1,H,N,Y,T,Z)

DIMENSION Y(N)

DOUBLE PRECISION X1,H,Y,T,Z,S,XX,P,Q

Z=0.0

IF (N.LE.0) RETURN

IF (N.EQ.1) THEN

Z=Y(1)

RETURN

END IF

IF (N.EQ.2) THEN

Z=(Y(2)*(T-X1)-Y(1)*(T-X1-H))/H

RETURN

END IF

IF (T.GT.X1) THEN

P=(T-X1)/H

I=P

Q=1.0*I

I=I+1

IF (P.GT.Q) I=I+1

ELSE

I=1

END IF

K=I-4

IF (K.LT.1) K=1

M=I+3

IF (M.GT.N) M=N

Z=0.0

DO 30 I=K,M

S=1.0

XX=X1+(I-1)*H

DO 20 J=K,M

IF (J.NE.I) THEN

XJ=X1+(J-1)*H

S=S*(T-XJ)/(XX-XJ)

END IF

20 CONTINUE

Z=Z+S*Y(I)

30 CONTINUE

RETURN

END

这是测试用的主程序代码。函数为F(x) = e^(-x)，求x = 0.63处的函数近似值：

PROGRAM EELGRO

DIMENSION Y(10)

DOUBLE PRECISION Y,H,X1,T,Z

DATA Y/0.904837,0.818731,0.740818,0.670320,0.606531,

* 0.548812,0.496585,0.449329,0.406570,0.367879/

H=0.1

N=10

X1=0.1

T=0.63

CALL EELGR(X1,H,N,Y,T,Z)

WRITE(*,20) T,Z

20 FORMAT(1X,'T=',F6.3,10X,'Z=',D15.6)

END

---