私のラグランジュ補間アルゴリズムが機能しないのはなぜですか？

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 

どこで失敗しましたか？

Lagrange interpolation method

事前にどうもありがとう：私は基本的にこれを実装する必要があります。

Answer 1

これは何をしますか？（完全な）プログラム

 real x = 1. 
     end 

を考えてみましょうか

これがフリーフォームソースの場合、無効なプログラムです。固定形式のソースであれば、それは有効なプログラムです。

固定形式のソースでは、列6の後のスペースはほとんど効果がありません。上記のプログラムは正確に

 realx=1. 
     end 

のようなもので、私たちはちょうど値1.を持つようにrealxと呼ばれる暗黙的に宣言した本当の変数を設定していることがわかります。

 implicit none 
     real x = 1. 
     end 

が表示されます。自由形式のソースで

、宣言文で初期化はそうのように、::を必要とします。

 real :: x = 1. 
     end 

そして：implicit noneを使用します。

Answer 2

（別の答えで説明されている）「シンボル連結」問題に加えて、は、I（各グリッドポイントのラグランジュ多項式を計算するため）ごとに1.0にリセットする必要があります。 （Y(I)）、その点の関数値で乗算します。

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 

：これら

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 

を固定した後、我々はおおよそ（=補間）と正確な結果とかなり良い合意を得ます