fortran - polinomio - método de interpolacion de lagrange

¿Por qué mi algoritmo de interpolación de Lagrange no funciona?

Question

Por alguna razón, nunca interpola, pero da 0 como respuesta. El código es:

  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

¿Dónde fallé? Básicamente necesito implementar esto:

Método de interpolación de Lagrange

Muchas gracias por adelantado.

Answer 1

Además del problema de "concatenación de símbolos" (explicado en la otra respuesta), parece que TEMP debe restablecerse a 1.0 por cada I (para calcular el polinomio de Lagrange para cada punto de la grilla), además tenemos que multiplicarlo por el valor funcional en ese punto ( Y(I) ). Después de arreglar estos

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

obtenemos una muy buena concordancia entre los resultados aproximados (= interpolados) y exactos:

 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 

Answer 2

Considere el programa (completo)

      real x = 1.
      end

¿Qué hace esto?

Si se trata de una fuente de forma libre, entonces es un programa no válido. Si es fuente de forma fija, entonces es un programa válido.

En fuente de forma fija, los espacios después de la columna 6 en gran medida no tienen ningún efecto. El programa de arriba es exactamente como

      realx=1.
      end

y podemos ver que solo estamos estableciendo una variable real implícitamente declarada llamada realx para tener el valor 1. .

      implicit none
      real x = 1.
      end

mostrará un problema

En la fuente de forma libre, la inicialización en una declaración de declaración requiere :: , así:

      real :: x = 1.
      end

Y: usar implicit none .