Why is my Lagrange interpolation algorithm not working?

Question

For some reason it never interpolates, but it gives 0 as an answer. The code is:

  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

Where did I fail? I basically need to implement this:

Lagrange interpolation method

Thanks a lot in advance.

Always use implicit none! Doubly so when having a problem and quadruply so when presenting the code to others but use it for the sanity of your own mind. — Vladimir F Героям слава
– Vladimir F Героям слава, Commented Oct 6, 2023 at 5:01

roygvib · Accepted Answer · 2017-05-04 18:29:25Z

In addition to the "symbol-concatenation" problem (explained in the other answer), it seems that TEMP needs to be reset to 1.0 for every I (to calculate the Lagrange polynomial for each grid point), plus we need to multiply it by the functional value on that point (Y(I)). After fixing these

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

we get a pretty good agreement between the approximate (= interpolated) and exact results:

 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

francescalus · Accepted Answer · 2023-10-05 23:17:33Z

2

Consider the (complete) program

      real x = 1.
      end

What does this do?

If this is free-form source then it is an invalid program. If it is fixed-form source then it is a valid program.

In fixed-form source, spaces after column 6 largely have no effect. The program above is exactly like

      realx=1.
      end

and we can see that we're just setting an implicitly declared real variable called realx to have value 1.. Such a statement is an assignment statement, not a type declaration statement with initialization.

Preventing implicit typing, as in

      implicit none
      real x = 1.
      end

will show a problem.

In both free- and fixed-form source, initialization in a declaration statement requires ::, like so:

      real :: x = 1.
      end

And: use implicit none.

edited Oct 5, 2023 at 23:17

answered May 4, 2017 at 16:50

francescalus

33k29 gold badges87 silver badges128 bronze badges

Collectives™ on Stack Overflow

Why is my Lagrange interpolation algorithm not working?

2 Answers 2

Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related