可查看visual fortran 常用算法程序集相關內容
拉格朗日插值
LAGINT.F
!
!根據給定點的函數值，應用拉格朗日插值公式，計算指定插值點處的
！函數值。
！在指定點的前后各取四個結點，用七次拉格朗日插值公式。
！X，Y為輸入參數，各有N個點
！T為指定插入點
！Z為輸出參數，是插值點T處的近似值。
SUBROUTIONE LAGINT（X，Y，N ，T，Z）
DIMENSION X（N），Y（N）
DOUBLE PRECISION X，Y，T，X，S
Z=0.0
IF (N.LE.0) RETURN
IF(N.EQ.1) THEN
Z=Y(1)
ENDIF
IF(N.EQ.2) THEN
Z=(Y(1)*(T-X(2))-Y(2)*(T-X(1)))/(X(1)-X(2))
RETURN
ENDIF
I=1
10 IF(X(1) .LT. T)
I=I+1
IF(I .LE. N) GOTO 10
ENDIF
K=I-4
IF(K .LT. 1) K=1
M=I+3
IF( M .GT. N) M=N
DO 30 I=K,M
S=1.0
DO 20 J=K,M
IF(J .NE. I) THEN
S=S*(T-X(J))/(X(I)-X(J))
ENDIF
20 CONTINUE
Z=Z+S*Y(I)
30 CONTINUE
RETURN
END
1.2 樣條插值
FORTRAN程序
SPLINE.F
! SPLINE.FOR -- Natural Cubic Spline
! Call SplineCalculate with the two arrays of control points (and
! the number of elements in each array). To interpolate the spline,
! Call SplineEvaluate for each 'X' for which you want a 'Y' along the 
! spline's curve.
! 

subroutine SplineCalculate( inx, ina, numx )
!ms$if .not. defined(LINKDIRECT)
!ms$attributes dllexport :: SplineCalculate
!ms$endif

implicit none
real*4 ina(*), inx(*)
integer numx

integer MAXPOINTS
parameter( MAXPOINTS = 1000 )
integer nxs
real*4 x(0:MAXPOINTS)
real*4 a(0:MAXPOINTS), b(0:MAXPOINTS)
real*4 c(0:MAXPOINTS), d(0:MAXPOINTS)
common /SplineData/ a, b, c, d, nxs, x

integer i
real*4 alpha(0:MAXPOINTS), l(0:MAXPOINTS)
real*4 mu(0:MAXPOINTS), z(0:MAXPOINTS)
real*4 h(0:MAXPOINTS)

! ===== Step 1
nxs = numx-1 
do i = 0, nxs
x(i) = inx(i+1)
a(i) = ina(i+1)
end do
do i = 0, nxs-1
h(i) = x(i+1) - x(i)
end do
x(nxs+1) = 40000

! ===== Step 2
do i = 1, nxs-1
alpha(i) = 3.0 * (a(i+1) * h(i-1)-a(i) * 
& (x(i+1)-x(i-1)) + a(i-1) * h(i)) / 
& (h(i-1) * h(i))
end do

! ===== Step 3
l(0) = 1.0
mu(0) = 0.0
z(0) = 0.0

! ===== Step 4
do i = 1, nxs-1
l(i) = 2.0 * (x(i+1)-x(i-1))-h(i-1) * mu(i-1)
mu(i) = h(i) / l(i)
z(i) = (alpha(i)-h(i-1) * z(i-1)) / l(i) 
end do

! ===== Step 5
l(nxs) = 1.0
z(nxs) = 0.0
c(nxs) = 0.0

! ===== Step 6
do i = nxs-1, 0, -1
c(i) = z(i) - mu(i)*c(i+1)
b(i) = (a(i+1)-a(i))/h(i) - h(i)* 
& (c(i+1)+2.0 * c(i)) / 3.0
d(i) = (c(i+1)-c(i)) / (3.0 * h(i)) 
end do

end

subroutine SplineEvaluate( evalx, evaly ) 
!ms$if .not. defined(LINKDIRECT)
!ms$attributes dllexport :: SplineEvaluate
!ms$endif

real*4 evalx, evaly

integer MAXPOINTS
parameter( MAXPOINTS = 1000 )
integer nxs
real*4 x(0:MAXPOINTS)
real*4 a(0:MAXPOINTS), b(0:MAXPOINTS)
real*4 c(0:MAXPOINTS), d(0:MAXPOINTS)
common /SplineData/ a, b, c, d, nxs, x

real*4 term
integer i

i = 0 
do while( .not. ((x(i) .le. evalx) .and. (x(i+1) .gt. evalx)) )
i = i + 1
end do 

term = evalx - x(i)
evaly = a(i) + b(i)*term + c(i)*term*term + 
& d(i)*term*term*term

end
1.3 曲線擬合
1.3.1.1 FORTRAN程序
!HPINT.FOR
!用最小二乘法求N個數據點的擬合多項式。
！X，Y 輸入參數，分別存放N個數據點的X坐標與Y坐標 
！A 輸出參數，存放M-1次擬合多項式的M個系數
！N 輸入參數，數據點數
！M 輸入參數，擬合多項式的項數，M N且M 20
！DT1，DT2，DT3—輸出參數， 分別為擬合多項式與數據點偏差的平方和、絕對
！ 值之和、絕對值的最大值
！
SUBROTINE HPINT(X,Y,A,N,M,DT1,DT2,DT3)
DIMENSION X(N), Y(N), A(M), S(20), T(20), B(20)
DOUBLE PRECISION X, Y, A, S, T, B, DT1, DT2, DT3, Z, D1, P,
& C, D2, G, Q, DT
DO 5 I=1,M
5 A(I)=0.0
IF ( M .GT. N) M=N
IF( M .GT. 20)M=20
Z=0.0
DO 10 I=1,N
10 Z=Z+X(I)/N
B(1)=1.0
D1=N
P=0.0
C=0.0
DO 20 I=1,N
P= P + (X(I) – Z)
C = C + Y(I)
20 CONTINUE
C = C/D1
P = P/D1
A(1) = C * B(1)
IF( M .GT. 1) THEN
T(2)=1.0
T(1) = -P
D2 = 0.0
C=0.0
G=0.0
DO 30 I=1,N
Q= X(I) –Z –P
D2=D2 + Q*Q
C= Y(I)*Q + C
G = (X(I) –Z)*Q*Q+G
30 CONTINUE
C = C/D2
P = G/D2
Q=D2/D1
D1=D2
A(2)=C*T(2)
A(1) = C*T(1) +A(1)
ENDIF
DO 100 J=3,M
S(J) = T(J-1)
S(J-1) = -P*T(J-1) + T(J-2)
IF ( J .GE. 4) THEN
DO 40 K=J-2,2,-1
40 S(K) = -P*T(K) + T(K-1) – Q*B(K)
ENDIF
S(1) = –P*T(1) – Q*B(1)
D2=0.0
C=0.0
G=0.0
DO 70 I=1,N
Q=S(J)
DO 60 K=J-1,1,-1
60 Q=Q* (x(I)-z_ + s(k)
D2= D2 + Q*Q
C= Y(I)*Q + C
G= ( X(I)-Z) QQ +G
70 CONTINUE
C C/D2
P=G/D2
Q=D2/D1
D1=D2
A(J) = C*S(J)
T(J) = S(J)
DO 80 K = J-1,1,-1
A(K) = C * S(K) +A(K)
B(K) = T(K)
T(K) = S(K)
80 CONTINUE
100 CONTINUE
DT1=0.0
DT2=0.0
DT3=0.0
DO 120 I=1,N
Q=A(M)
DO 110 K=M-1,1,-1
110 Q= Q*(X(I)-Z) + A(K)
DT=Q-Y(I)
IF (ABS(DT) .GT. DT3) D3=ABS(DT)
DT1= DT1+DT*DT
DT2=DT2+ABS(DT)
120 CONTINUE
RETURN
END