[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Help-gsl] Incomplete elliptic integral E(phi, k) decreasing with ph
From: |
Liam Healy |
Subject: |
Re: [Help-gsl] Incomplete elliptic integral E(phi, k) decreasing with phi? |
Date: |
Sat, 22 Mar 2008 17:15:55 -0400 |
Frank,
I have written a little C program to show this:
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_mode.h>
#include <gsl/gsl_sf_ellint.h>
int
main (void)
{
double phi = 0.6*M_PI;
double k = 0.5;
gsl_sf_result result;
int status = gsl_sf_ellint_E_e (phi, k, GSL_PREC_DOUBLE, &result);
printf ("status = %s\n", gsl_strerror(status));
printf ("answer = %.18f\n"
" +/- % .18f\n",
result.val, result.err);
return status;
}
I compiled it with
gcc specfun_e.c -lgsl -lgslcblas -o ellint
When run, I get
status = success
answer = 1.193936630679625743
+/- 0.000000000000000647
When I change the value of phi to M_PI, I get
status = success
answer = 0.000000000000000122
+/- 0.000000000000000000
I'm not sure what you mean by "seem to be correct also";
do you mean you got my results, or you got the correct results?
Another set of values which seem wrong:
phi = -0.5*M_PI;
gives
answer = 1.467462209339427393
while
phi = -0.49999999*M_PI;
gives
answer = -1.467462183529858910
The positive number is wrong; I am sure that this function is
continuous in any case and shouldn't show a jump at -pi/2.
By the way I am using GSL 1.8 in Debian stable; if you get the correct
answers is it possible something was fixed in subsequent versions?
Thanks,
Liam
On 3/22/08, Frank Reininghaus <address@hidden> wrote:
> Hi Liam,
>
>
> On Saturday 22 March 2008 20:15:06 Liam Healy wrote:
> > Yet I try gsl_sf_ellint_E_e
> > for phi = 0.5pi and k=0.5, I find the value returned is 1.46746 (which
> > agrees with the result from the complete elliptic integral as it
> > should),
> > and for phi= 0.6pi and k=0.5, value is 1.19394. In fact, for phi=pi,
> > I get essentially 0, when it looks like I should get 2*1.46746 because
> > sin^2 is symmetric about pi/2.
>
>
> I tried what you decribed and got results that seem to be correct also for
> phi=0.6*pi and phi=pi. Could you send a complete (but short) program that
> shows the behavior you mentioned?
>
> Thanks,
>
> Frank
>
>
>