Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Error writing recursive function for Bonnet's recursion formula in C

  1. Jan 9, 2012 #1
    Bonnet’s recursion formula for Legendre polynomials is
    P(n,x)=1 for n=0
    P(n,x)=x for n=1 and
    P(n,x)= (2n-1)/n*x*P(n-1,x)-(n-1)/n*P(n-2,x)

    I tried to write recursive function in C to calculate the value of polynomial for a given n and x
    Code (Text):
    float legpol(int n, float x)
        {
            if(n==0)
                return 1;
            else if(n==1)
                return x;
            else
                return(2*n-1)/n*x*legpol(n-1,x)-(n-1)/n*legpol(n-2,x);
        }
     
    But this doesn't give the correct result. What's wrong with the function??

    Thanks in advance
     
  2. jcsd
  3. Jan 9, 2012 #2
    Can you give a sample of the output vs. what the correct answer should be?
     
  4. Jan 9, 2012 #3
    The actual legendre's polynomials are
    for n = 0 1
    for n= 1 x
    for n = 2
    82722ebdf5cc1c9deb3475a7ec25bd3d.png
    for n= 3
    b8e8e3c74d685b1d03f1d11dd2b0a78d.png
    for n =4
    1b43b5affe459365b174311803a08ad6.png
    and so on


    if n= 2 and x= .1 the output must be -0.485 but my recursive function gives 0.01 which is wrong
     
  5. Jan 9, 2012 #4

    Mark44

    Staff: Mentor

    The problem is, I believe, that your formula is using int division in a couple of places, and they give you incorrect function values. Inexperienced C and C++ programmers often don't realize that there are two kinds of division - integer and floating point.

    Your code has this line:
    Code (Text):
    return(2*n-1)/n*x*legpol(n-1,x)-(n-1)/n*legpol(n-2,x);
    If n == 3, for example, the (2*n - 1)/n expression evaluates to 5/3 == 1. The other int term is multiplied by 2/3 == 0.

    An easy fix is to change the line above to
    Code (Text):
    return(2*n-1.0)/n*x*legpol(n-1,x)-(n-1.0)/n*legpol(n-2,x);
    This change causes the numerators of both fractions to be promoted to double before the division actually occurs, forcing floating point division to be performed, rather than integer division.
     
  6. Jan 9, 2012 #5
    Thank you!! This change worked fine.

    I tried to figure out the error and checked the value of 'n' after each function call and thought there might have been changes in the value of 'n' within the recursive function call. So i tried to preserve the value of 'n' by creating another variable 'b' just to use it in nested function call. I modified the above function like this::
    Code (Text):

    float legpol(int n, float x)
        {
            int b;
            if(n==0)
                return 1;
            else if(n==1)
                return x;
            else
            {
                 b=n;
                 return(2*b-1)/b*x*legpol(b-1,x)-(b-1)/b*legpol(b-2,x);
             }
         }
     
    Amazingly it worked. In this case why shouldn't i write (2*b-1.0) ??
     
  7. Jan 9, 2012 #6

    Mark44

    Staff: Mentor

    I don't know why your revision should work. For some int divisions you get the answer you would expect, such as 6/3 or 10/2. It might be that what got evaluated was something like this.

    In any case, I would get rid of the b variable and use code that forces floating point division, such as 2nd example in my previous post.
     
  8. Jan 10, 2012 #7
    I'm sorry for the second revision i posted. I had actually typed
    Code (Text):
    float b;
    in the third line of my function not
    Code (Text):
     int b
    .

    I now know what the problem was. Thank you very much for the solution. It really helped me.
     
  9. Jan 10, 2012 #8

    Mark44

    Staff: Mentor

    Yes, that would have made a difference, by causing floating point division to be performed instead of integer division.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Error writing recursive function for Bonnet's recursion formula in C
  1. Recursion in C (Replies: 12)

Loading...