Res: Res: [Help-gsl] Announcement and Question - jacobi
Thanks! I will take a closer look and will implement this in the following days.
Paulo
- Mensagem original
De: Brian Gough <[EMAIL PROTECTED]>
Para: help-gsl@gnu.org
Cc: Paulo Jabardo <[EMAIL PROTECTED]>
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.c3 Jul 2005 10:29:26 -1.37
+++ legendre_poly.c5 Dec 2006 19:50:09 -
@@ -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
___
Help-gsl mailing list
Help-gsl@gnu.org
http://lists.gnu.org/mailman/listinfo/help-gsl
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 - 1.37
+++ legendre_poly.c 5 Dec 2006 19:50:09 -
@@ -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
Res: [Help-gsl] Announcement and Question - jacobi
That's exactly the one I tried to "reverse engineer" since Jacobi polynomial are small modification of that. Paulo - Mensagem original De: Brian Gough <[EMAIL PROTECTED]> Para: help-gsl@gnu.org Cc: Paulo Jabardo <[EMAIL PROTECTED]> Enviadas: Segunda-feira, 4 de Dezembro de 2006 14:53:03 Assunto: Re: [Help-gsl] Announcement and Question - jacobi At Fri, 1 Dec 2006 03:38:06 -0800 (PST), Paulo Jabardo wrote: > I still have not implemented the error estimates. I tried to look at > the code used to compute Legendre polynomials but I couldn't figure > where did the error come from. Apparently it is related to the > number of floating point operations. Does any one have any hints or > references on that that I could use? Which function were you looking at - gsl_sf_legendre_Pl_e? -- Brian Gough ___ Novidade no Yahoo! Mail: receba alertas de novas mensagens no seu celular. Registre seu aparelho agora! http://br.mobile.yahoo.com/mailalertas/ ___ Help-gsl mailing list Help-gsl@gnu.org http://lists.gnu.org/mailman/listinfo/help-gsl
Re: [Help-gsl] Announcement and Question - jacobi
At Fri, 1 Dec 2006 03:38:06 -0800 (PST), Paulo Jabardo wrote: > I still have not implemented the error estimates. I tried to look at > the code used to compute Legendre polynomials but I couldn't figure > where did the error come from. Apparently it is related to the > number of floating point operations. Does any one have any hints or > references on that that I could use? Which function were you looking at - gsl_sf_legendre_Pl_e? -- Brian Gough ___ Help-gsl mailing list Help-gsl@gnu.org http://lists.gnu.org/mailman/listinfo/help-gsl
[Help-gsl] Announcement and Question - jacobi
Hello, I just released an extension to gsl that computes jacobi polynomials, zeros of jacobi polynomials and several operations related to Gauss-Jacobi quadrature (numerical integration, numerical derivatives and interpolation using the quadrature nodes). The extension can be found at: http://strouhal.mcef.ep.usp.br/pg/pjabardo/jacobi Now comes the question: I still have not implemented the error estimates. I tried to look at the code used to compute Legendre polynomials but I couldn't figure where did the error come from. Apparently it is related to the number of floating point operations. Does any one have any hints or references on that that I could use? Thanks for any help. Paulo Jabardo ___ Novidade no Yahoo! Mail: receba alertas de novas mensagens no seu celular. Registre seu aparelho agora! http://br.mobile.yahoo.com/mailalertas/ ___ Help-gsl mailing list Help-gsl@gnu.org http://lists.gnu.org/mailman/listinfo/help-gsl
