FORTRAN - Neville's algorithm
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
INTEGER j,k
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 DO
END DO
res=P(1)
err=P(1)-P(2)
return
END SUBROUTINE Neville
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
INTEGER j,k
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 DO
END DO
res=P(1)
err=P(1)-P(2)
return
END SUBROUTINE Neville
ASKER CERTIFIED SOLUTION
membership
This solution is only available to members.
To access this solution, you must be a member of Experts Exchange.
ASKER
Numerical Recipes is actually a book from where I study numerical methods..
I note that you are starting your counters at 1. According to the following link, the algorithm was intended to be started from index 0, and suggests that if one wishes to start from 1, then of course adjustment to the boundary conditions at the top needs to be made.
http://www.cs.uaf.edu/~bueler/nevM.htm
If that's not the issue (I've only glanced at the coding), then let us know
http://www.cs.uaf.edu/~bueler/nevM.htm
If that's not the issue (I've only glanced at the coding), then let us know
SOLUTION
membership
This solution is only available to members.
To access this solution, you must be a member of Experts Exchange.
ASKER
Basicly with my Subroutine I get the same results as with the one written in Numerical Recipes F77, but I0m curious can someone test my subroutine with the function 1/sin(x) and others ... I mean write a main program and tell me their experience?... Have I really done something wrong or do the results really look that odd as mines.
moorhouselondon:
Unlike C or C++, Fortran starts it's counters(indexes) for Arrays with number 1 not 0. for example Array[1] is actually the first member of the array, not the second like in C or C++ .. like I said
moorhouselondon:
Unlike C or C++, Fortran starts it's counters(indexes) for Arrays with number 1 not 0. for example Array[1] is actually the first member of the array, not the second like in C or C++ .. like I said
ASKER
Figuered thatone out ... the main problem was the function that I was trying to interpolate ... soo the algorithm is basicly ok! But the function is not! :((
Thanx for trying guys!
Thanx for trying guys!
ASKER