为什么我的拉格朗日插值算法不起作用？

Question

由于某种原因，它从不插值，但它给出 0 作为答案。代码是：

  PROGRAM LAGRANGE
    REAL X(0:100), Y(0:100), INTERP
    REAL TEMP = 1.0
    REAL POLINOM = 0.0
    N=10
    OPEN(1,FILE="datos.txt")
    DO I=0,100 !We 'clean' the arrays: all positions are 0
       X(I)=0.0
       Y(I)=0.0
    END DO        
    DO I=0,10 !We read the data file and we save the info
       READ(1,*) X(I), Y(I)
    END DO
    CLOSE(1)

    WRITE(*,*) "Data table:"
    DO I=0,10
    WRITE(*,*) X(I), Y(I)
    END DO


    WRITE(*,*) "Which value of X do you want to interpolate?"
    READ(*,*) INTERP

    DO I=0,N
        DO J=0,N
            IF(J.NE.I) THEN !Condition: J and I can't be equal
                TEMP=TEMP*(INTERP-X(J))/(X(I)-X(J))
            ELSE IF(J==I) THEN
                TEMP=TEMP*1.0
            ELSE
            END IF
        END DO
        POLINOM=POLINOM+TEMP
    END DO

    WRITE(*,*) "Value: ",POLINOM

  STOP
  END PROGRAM

我哪里失败了？我基本上需要实现这个：

提前非常感谢。

Answer 1

除了“符号串联”问题（在另一个答案中解释）之外，似乎每个

TEMP

都需要重置为 1.0 （以计算每个网格点的拉格朗日多项式），加上我们需要将其乘以该点的函数值 (

I

)。修复这些之后Y(I)

我们在近似值（=插值）和精确结果之间得到了很好的一致性：

PROGRAM LAGRANGE
    implicit none  !<-- always recommended
    REAL :: X(0:100), Y(0:100), INTERP, TEMP, POLINOM
    integer :: I, J, K, N

    N = 10
    X = 0.0
    Y = 0.0

    !! Test data (sin(x) over [0,2*pi]).
    DO I = 0, N
        X(I) = real(I) / real(N) * 3.14159 * 2.0
        Y(I) = sin( X(I) )
    END DO

    WRITE(*,*) "Data table:"
    DO I = 0, N
        WRITE(*,*) X(I), Y(I)
    END DO

    interp = 0.5   !! test value

    POLINOM = 0.0
    DO I = 0, N

        TEMP = 1.0   !<-- TEMP should be reset to 1.0 for every I
        DO J = 0, N
            IF( J /= I ) THEN
                TEMP = TEMP * (interp - X(J)) / (X(I) - X(J))
            END IF
        END DO
        TEMP = TEMP * Y(I)  !<-- also needs this

        POLINOM = POLINOM + TEMP
    END DO

    print *, "approx : ", POLINOM
    print *, "exact  : ", sin( interp )
end

Answer 2

 Data table:
   0.00000000       0.00000000    
  0.628318012      0.587784827    
   1.25663602      0.951056182    
   1.88495409      0.951056957    
   2.51327205      0.587786913    
   3.14159012       2.53518169E-06
   3.76990819     -0.587782800    
   4.39822626     -0.951055467    
   5.02654409     -0.951057792    
   5.65486193     -0.587789178    
   6.28318024      -5.07036339E-06
 approx :   0.479412317    
 exact  :   0.479425550 

这个有什么作用？

如果这是自由格式的源代码，那么它是一个无效的程序。如果它是固定格式的源代码，那么它就是一个有效的程序。

在固定格式的源代码中，第 6 列之后的空格基本上没有影响。上面的程序一模一样

      real x = 1.
      end

我们可以看到，我们只是将一个名为 

      realx=1.
      end

的隐式声明实数变量设置为具有值

realx

。这样的语句是赋值语句，而不是带有初始化的类型声明语句。防止隐式键入，如

1.

会显示问题。

在自由格式和固定格式源代码中，声明语句中的初始化需要

      implicit none
      real x = 1.
      end

，如下所示：

::

并且：使用 

      real :: x = 1.
      end

。