Re: Res: [Help-gsl] Announcement and Question - jacobi

[Brian Gough]

At Mon, 4 Dec 2006 09:59:05 -0800 (PST),
Paulo Jabardo wrote:
> 
> That's exactly the one I tried to "reverse engineer" since Jacobi
> polynomial are small modification of that.
> 

For the recurrence relation, the approximation there is very crude
just O(number operations) * result * GSL_DBL_EPSILON.  That will not
be reliable when there is cancellation in the recurrence.

In general, to calculate the error reliably it needs to be propagated
with the computed value.  The example below shows the simple way to do
it.  In this case the error recurrence could be simplified by making
the approximation ell>>1.

-- 
Brian Gough

Network Theory Ltd,
Publishing Free Software Manuals --- http://www.network-theory.co.uk/

Index: legendre_poly.c
===================================================================
RCS file: /home/gsl-cvs/gsl/specfunc/legendre_poly.c,v
retrieving revision 1.37
diff -u -r1.37 legendre_poly.c
--- legendre_poly.c     3 Jul 2005 10:29:26 -0000       1.37
+++ legendre_poly.c     5 Dec 2006 19:50:09 -0000
@@ -141,16 +141,25 @@
     double p_ellm2 = 1.0;    /* P_0(x) */
     double p_ellm1 = x;      /* P_1(x) */
     double p_ell = p_ellm1;
+
+    double e_ellm2 = GSL_DBL_EPSILON;
+    double e_ellm1 = fabs(x)*GSL_DBL_EPSILON;
+    double e_ell = e_ellm1;
+
     int ell;
 
     for(ell=2; ell <= l; ell++){
       p_ell = (x*(2*ell-1)*p_ellm1 - (ell-1)*p_ellm2) / ell;
       p_ellm2 = p_ellm1;
       p_ellm1 = p_ell;
+
+      e_ell = fabs(x*(2*ell-1.0)/ell) * e_ellm1 + fabs((ell-1.0)/ell)*e_ellm2;
+      e_ellm2 = e_ellm1;
+      e_ellm1 = e_ell;
     }
 
     result->val = p_ell;
-    result->err = (0.5 * ell + 1.0) * GSL_DBL_EPSILON * fabs(p_ell);
+    result->err = e_ell;
     return GSL_SUCCESS;
   }
   else {