Advanced Search
All Posts : Product Tips & Techniques

Nice numerical trouble

Axel Vogt 5821 Maple
December 04 2006
1
Playing with numerics for the hypergeometric function 2F1 i stumbled into the following exactness stuff, where lots of digits are lost and one needs to increase working precision quite a lot to get a good answer: f1 and f2 result from evaluating ugly, but usual transcendent functions, f3 comes through high precision and f0 is the limiting case. 
  restart;
  
  hypergeom([1,2+epsilon],[3],z): 
  '%' = simplify(convert(%,StandardFunctions));
  theSol1:=rhs(%): 
  
  theSol1 = -2*(z*(1-z)^epsilon*epsilon-1+(1-z)^epsilon)/
    z^2/((1-z)^epsilon)/epsilon/(1+epsilon);
  is(%);
  theSol2:=rhs(%%):
  
  limit(hypergeom([1, 2+epsilon],[3],z),epsilon = 0):
  '%'= simplify(%);
  theSol0:=rhs(%):

  hypergeom([1, 2 + epsilon], [3], z) =

                     (-epsilon)
          2 (-(1 - z)           + 1 + epsilon z)
        - --------------------------------------
                  2
                 z  (1 + epsilon) epsilon


               (-epsilon)
    2 (-(1 - z)           + 1 + epsilon z)
  - -------------------------------------- =
            2
           z  (1 + epsilon) epsilon

                      epsilon                      epsilon
          2 (z (1 - z)        epsilon - 1 + (1 - z)       )
        - -------------------------------------------------
                2        epsilon
               z  (1 - z)        epsilon (1 + epsilon)


                                 true


      lim       hypergeom([1, 2 + epsilon], [3], z) =
  epsilon -> 0

          2 (z + ln(1 - z))
        - -----------------
                  2
                 z

  gc();
  remDigits:=Digits: Digits:=18;
  Epsilon:=10^(-15):
  tstData:=[a=1, b=2 + Epsilon, c=3,  z= 3/2];
  ``;
  subs(epsilon=Epsilon,theSol1): eval(%,tstData): f1:=evalf(%);
  
  subs(epsilon=Epsilon,theSol2): eval(%,tstData): f2:=evalf(%);
  
  Digits:=2*Digits:
  subs(epsilon=Epsilon,theSol1): eval(%,tstData): f:=evalf(%): 
  Digits:= floor(Digits/2): f3:=evalf(f);
  Digits:=remDigits:
  
  subs(epsilon=Epsilon,theSol0): eval(%,tstData): f0:=evalf(%);
  
  Digits:=remDigits:

                             Digits := 18


                              2000000000000001
       tstData := [a = 1, b = ----------------, c = 3, z = 3/2]
                              1000000000000000


                 f1 := -0.720 - 2.79252679411846980 I
         f2 := -0.711111111111119664 - 2.79252680319092647 I
         f3 := -0.717202506168940958 - 2.79252680319092647 I
         f0 := -0.717202506168937500 - 2.79252680319092732 I
convert

Please Wait...