Legendre polynomials in the reverse direction

[ognik]

I have just written a program to calculate Legendre Polynomials, finding for P_l+1 using the recursion (l+1)P_l+1 + lP_l-1 - (2l+1).x.P_l=0 That is working fine.
The next section of the problem is to investigate the recursive polynomial in the reverse direction. I would solve this for P_l-1 in terms of P_l and P_l+1 - could someone please confirm if that is the correct way to 'reverse'? Also I would use the Rodriques formula to get the 2 starting values of P_n and P_n-1, it's a complex formula, is there a better way to get those starting values? Thanks for all help.

 

[ognik]

Hi again - I found enough to convince me that solving for P_l-1 is the correct way to calculate in the reverse direction. I also decided against using the Rodriguez formula, since I work out P_l and P_l+1 in the forward program, I can just use those as starting values for the reverse direction. All seemed good. However what I expect in the reverse direction, is to get P(x) for l=1, which should be always x. I figured that the last 2 P values correspond to L_input and L_input - 1, so loop down from L_input - 2 down to l=1 - see program code below.
I have spent hours with the debugger and pouring over the code, can't spot the bug or mistake in approach. Help please?

program LegendreDual
   implicit none
   integer, parameter   ::   rep=60
   integer, parameter   ::   dp = selected_real_kind(15, 307)
   real(kind=dp)     ::   X,P,P_minus,P_plus,saveP
   integer         ::   j, L_input
   character(len=rep )   ::   stars=REPEAT('*',rep)
   !----- header ----   
   print *, " "
   print *, 'Program LegendreDual starts'
   print *, stars 
   !----- Initialise ----   
   print *, " "
   print *,'Enter  l, X (l .LT.  0 to stop)'
   read *, L_input, X
   print *, " "
  while (L_input.GE.0) do
       print 20, 'L input','X Input','P'
     If (L_input.EQ.0) then
     P=0
  else if (L_input.EQ.1) then
       P=X
  else
     P_minus=1
  P=X
  Do j=1, L_input-1,+1           ! loop for fwd recursion
     P_plus=(((2*j)+1)*X*P)-(j*P_minus)
     P_plus=P_plus/(j+1)
         saveP=P
     P_minus=P                 ! roll values
     P=P_plus
       end do
     Print 50, L_input, X, P   
  P=saveP                   ! P from fwd; we can use P_plus and P
!-----------------------               ! to start the backward recursion
      Do j=L_input-2,1,-1               ! loop for reverse recursion,
     P_minus=(((2*j)+1)*X*P)-((j+1)*P_plus)   ! using P(n+1) and P just found
     P_minus=P_minus/j
     P=P_minus                   ! roll values
     P_plus=P
       end do     
  end if
  Print 50, j, X, P             ! only print once if l=0 or 1
     print *, stars  
     print *, " "
     print *,'Enter  l, X (l .LT.  0 to stop)'
     read *, L_input, X
  end do
   print *, " "
!--- admin ----
   print *, 'User selected end'
   print *, stars  
20   format (3X, A, 3X, A, 18X,A)
50   format (3X, I3, 5X, F10.5, 5X, F25.5)
end program LegendreDual

 

[ognik]

Please feel free to close or delete this thread and sorry for any inconvenience