2017-05-04 34 views
3

由於某種原因,它永遠不會插值,但它將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

非常感謝提前。

回答

2

考慮(完整的)程序

 real x = 1. 
     end 

這是什麼呢?

如果這是自由格式的源,那麼它是一個無效的程序。如果它是固定源,那麼它是一個有效的程序。

在固定格式源中,第6列之後的空格很大程度上不起作用。上面的程序是完全一樣

 realx=1. 
     end 

,我們可以看到,我們只是設置所謂realx擁有價值1.隱式聲明的真正變量。

 implicit none 
     real x = 1. 
     end 

會出現問題。

在自由格式源,初始化的聲明語句需要::,就像這樣:

 real :: x = 1. 
     end 

而且:使用implicit none

3

除了「符號並置」的問題(在對方的回答解釋),似乎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 

我們得到近似的(=插值)和詳細的結果之間一個很好的協議:

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