help-gsl
[Top][All Lists]
Advanced

[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
>
>
>




reply via email to

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