help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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


From: Paulo Jabardo
Subject: Res: Res: [Help-gsl] Announcement and Question - jacobi
Date: Tue, 5 Dec 2006 13:04:05 -0800 (PST)

Thanks! I will take a closer look and will implement this in the following days.

Paulo


----- Mensagem original ----
De: Brian Gough <address@hidden>
Para: address@hidden
Cc: Paulo Jabardo <address@hidden>
Enviadas: Terça-feira, 5 de Dezembro de 2006 18:02:27
Assunto: Re: Res: [Help-gsl] Announcement and Question - jacobi

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 {






                
_______________________________________________________ 
O Yahoo! está de cara nova. Venha conferir! 
http://br.yahoo.com




reply via email to

[Prev in Thread] Current Thread [Next in Thread]