Why is my Lagrange interpolation algorithm not working?
For some reason, it never interpolates, but gives 0 as an answer. Code:
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:
Thank you very much in advance.
source to share
In addition to the "character concatenation" problem (explained in another answer), it seems like there TEMP
should be a reset to 1.0 for each I
(to calculate the Lagrange polynomial for each grid point), plus we need to multiply it by the function value at 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 pretty good agreement between approximate (= interpolated) and precise 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
source to share
Consider the (complete) program
real x = 1.
end
What does it do?
If it is a free-form source, it is an invalid program. If it is a fixed-form source, then it is a valid program.
In a fixed form source, the spaces after column 6 are largely unaffected. The program above exactly matches
realx=1.
end
and we see that we are just setting the implicitly declared real variable named realx
to have a value 1.
.
implicit none
real x = 1.
end
will show the problem.
In a free-form source, initialization in a declaration statement requires ::
, for example:
real :: x = 1.
end
And: use implicit none
.
source to share