FORTRAN - Neville's algorithm
Posted on 2004-11-07
Hello! Don't see a Fortran Language topic so I'm posting here! Soo
I'm Having a lil' problem here ... Trying to write a good Interpolation module using Neville's algorithm and it just seems like I can't get a good interpolation for the function 1/sin(x).
Anyone have any ideas how to make an alternative subroutine or maybe what's wrong with mine? I'm posting here mine.
I know that interpolation 1/sin(x) is tricky beacuse the function has breaks around 0, pi, and 2 pi, but nevertheles
the numbers I'm getting have no sense
I'm trying to interpolate with 10th order polinomial and trying to interpolate 30 equidistant points in interval (0, 2*pi)
and all I'm getting out are really wierd numbers ... am I doing something wrong?
Please Help ...
SUBROUTINE Neville(A, B, xm, xt, res, err)
REAL, dimension(:), intent( in ) :: A, B
REAL, intent( in ) :: xt
REAL, intent( out ) :: res, err
INTEGER, intent ( in ) :: xm
REAL, dimension(xm) :: P
P(:) = B(:)
DO j=1, xm
DO k=1, xm - j
P(k)=((xt - A(k+j)) * P(k) + (A(k)-xt) * P(k+1))/(A(k)-A(k+j))
END SUBROUTINE Neville