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 { _______________________________________________ Help-gsl mailing list Help-gsl@gnu.org http://lists.gnu.org/mailman/listinfo/help-gsl
