FORTRAN - Neville's algorithm

Posted on 2004-11-07
Last Modified: 2013-11-08
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

Question by:Physicistm
    LVL 3

    Author Comment

    Ohh yes by the way ... It's Fortran 90!!
    LVL 55

    Accepted Solution

    Have a look to this page:
    LVL 3

    Author Comment

    Numerical Recipes is actually a book from where I study numerical methods..
    LVL 31

    Expert Comment

    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.

    If that's not the issue (I've only glanced at the coding), then let us know
    LVL 31

    Assisted Solution

    The link cited by jaime_olivares uses 1-based indexing, so I think my initial thought is wrong.
    LVL 3

    Author Comment

    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.

    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
    LVL 3

    Author Comment

    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!

    Featured Post

    Free Trending Threat Insights Every Day

    Enhance your security with threat intelligence from the web. Get trending threat insights on hackers, exploits, and suspicious IP addresses delivered to your inbox with our free Cyber Daily.

    Join & Write a Comment

    How to remove superseded packages in windows w60 or w61 installation media (.wim) or online system to prevent unnecessary space. w60 means Windows Vista or Windows Server 2008. w61 means Windows 7 or Windows Server 2008 R2. There are various …
    The purpose of this article is to demonstrate how we can use conditional statements using Python.
    This theoretical tutorial explains exceptions, reasons for exceptions, different categories of exception and exception hierarchy.
    In this seventh video of the Xpdf series, we discuss and demonstrate the PDFfonts utility, which lists all the fonts used in a PDF file. It does this via a command line interface, making it suitable for use in programs, scripts, batch files — any pl…

    755 members asked questions and received personalized solutions in the past 7 days.

    Join the community of 500,000 technology professionals and ask your questions.

    Join & Ask a Question

    Need Help in Real-Time?

    Connect with top rated Experts

    18 Experts available now in Live!

    Get 1:1 Help Now