如何用fortran编写一个线性插值程序

这是一元全区间等距插值子程序，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

去年学的Fortran语言，许多算法都忘仔野竖了。给你一个去年写的牛顿插值程序念大，需要改动的地方自己改，貌似循环计算（多层计算）那块处理得不错。

real*8 x(0:5),y(0:5),c(0:5),u1,u2,ci

data y/0.796,0.773,0.744,0.704,0.656,0.602/

open(1,file='7-2.dat')

write(*,*)'input u1=?,u2=?'

read(*,*)u1,u2

do i=0,5

x(i)=(i+1)*0.125

enddo

c(0)=y(0)

do i=1,5

call chashang(x,5,i,ci)

c(i)=ci

write(*,*)'c(i)=,i=',c(i),i

write(1,*)'脊羡c(i)=,i=',c(i),i

enddo

call newchazhi(u1,c,x,5)

call newchazhi(u2,c,x,5)

end

subroutine newchazhi(u,c,x,m)

real*8 c(0:m),x(0:m),u,N

S=c(0)

do i=1,5

N=c(i)

do j=0,i-1

N=N*(u-x(j))

enddo

S=S+N

enddo

write(*,*)u,s

write(1,*)u,s

end

subroutine chashang(x,m,i,ci)

real*8 y(0:5),x(0:m),f(1:i),ci

y(0)=0.796

y(1)=0.773

y(2)=0.744

y(3)=0.704

y(4)=0.656

y(5)=0.602

k=0

10 do j=1,i-k

f(j)=(y(j)-y(j-1))/(x(j+k)-x(j-1))

enddo

if(j==1) then

goto 5

else

do j=1,i-k

y(j-1)=f(j)

enddo

k=k+1

goto 10

endif

5 ci=f(j)

return

end

f95的侍散数。

彭的书上现成的程序。包括用sgl画图部分，函数老首可以自己改成多项式形式的，可以参考一下：

module INTERPOLATE_UTILITY

use sgl

implicit none

type point

real x,y

end type

real, parameter :: PI=3.14159

real, parameter :: xmin = 0.0, xmax = PI*3.0

integer, parameter :: N = 10, NP = 100

type(point), save :: datas(N)

type(point), save :: interpolate(NP)

real, save :: table(N,N), width

contains

! 产生数列

subroutine GenerateData(func)

real, external :: func

real r

integer i

width = (xmax-xmin)/(N-1)

r = 0

do i=1,N

datas(i)%x = r

datas(i)%y = func(r)

r = r+width

end do

end subroutine

! 建立difference table

subroutine BuildTable()

integer row,col,i

!real table(N,N)

table = 0

do i=1,N

table(i,1) = datas(i)%y

end do

do col=2,N

do row=1,N-col+1

table(row,col) = table(row+1, col-1) - table(row, col-1)

end do

end do

end subroutine

real function newton(x, th, num)

real x

integer th, num

real s, sum, coeff

integer f,i,j

if ( th+num-1 >N ) then

write(*,*) "数据点不足"

return

end if

newton = table(th,1)

s = (x-datas(th)%x)/width

f = 1

coeff = 1.0

do i=1,num-1

f = f*i

coeff = coeff*(s-i+1)

newton = newton + coeff*table(th,i+1)/real(f)

end do

end function

! 绘图函数

subroutine display()

real, parameter :: size = 0.1

integer i

call sglClearBuffer()

call sglColor3i(255,255,255)

! 把所有插值出来的点用线段掘历连接起来

do i=1,NP-1

call sglLineV( interpolate(i)%x, interpolate(i)%y,&

interpolate(i+1)%x, interpolate(i+1)%y)

end do

call sglColor3i(255,0,0)

! 画出n个数据点的位置

do i=1,N

call sglLineV( datas(i)%x-size, datas(i)%y-size,&

datas(i)%x+size, datas(i)%y+size)

call sglLineV( datas(i)%x+size, datas(i)%y-size,&

datas(i)%x-size, datas(i)%y+size)

end do

call sglUpdateBuffer()

end subroutine

end module

program main

use INTERPOLATE_UTILITY

implicit none

real, intrinsic :: sin

real xinc,x

integer i

call GenerateData(sin) ! 产生数据点

call BuildTable()

x=0

xinc = (xmax-xmin)/(NP-1)

do i=1,NP

interpolate(i)%x = x

interpolate(i)%y = newton(x,1,N) ! 插值出f(x)的值

x = x+xinc

end do

call sglDisplaySub(display)

call sglSetVirtual(xmin, 2.0, xmax, -2.0)

call sglCreateWindow(100,100,400,400,1)

call sglMainLoop()

stop

end program

---