Home › Forums › Student Help Forum › Math Advice (Student)

Subtraction of two generalized hypergeometric functions

Posted on 2008-07-17 04:28 By vishwa_r (8)
Categories:

Greetings to all,

I am trying to evaluate anexpression of the form :   P * xa * 1F2(a; 1+a-b; 1+a; x)  -  Q * xb  * 1F2(b; 1+b-a; 1+b; x)   where P,Q,a and b are constants. I need to evalute this expression for large values of the parameter x ( in loop). Eventhough the value of the above expression is a small quantity nearing unity, individual values of the hypergeometric functions will be very huge, and hence can create problems with the resolution of the processing. One way is to change the resolution to some huge value on case by case basis. Is there any trick or method so that it may be possible to evaluate this expression using the smaller fixed resolution values?. Thank you very much for any help.

Best regards,

Vishwa 

 

 

‹ Help! Is this right? statistics ›
alec's picture

Asymptotic

Comment on 2008-07-17 05:13 from alec (Maple Rating 4 1443)

It looks as if the 2 functions have different asymptotics, so their values shouldn't be very close for large x with constant P, Q, a, and b. Do you have a particular example?

Alec

Expression with typical values

Comment on 2008-07-17 06:40 from vishwa_r (8)

Hi alec,

 Thanks for the reply. I have attached the maple worksheet with the details and the reason for the values to be expected to be closer. Please have a look. I was wondering is there anyway to get rid of the huge variable resolution issue. Thanx for the help.

Best regards,

Vishwa

View 8014_Hypergeometric-diff.mw on MapleNet or Download 8014_Hypergeometric-diff.mw
View file details

 

acer's picture

possibility

Comment on 2008-07-17 10:05 from acer (Maple Rating 4 1106)

I figure that by "resolution" you mean Digits, the setting of the precision for floating-point calculation.

I'm not sure whether the following is useful to you,


> restart:
>
> p := Pi*csc(Pi*alp)
> *((E*x)^m*hypergeom([m],[1-alp,1+m],E*x)
> /(GAMMA(c)*GAMMA(1-alp)*GAMMA(1+m))
> -(E*x)^c*hypergeom([c], [1+alp, 1+c], E*x)
> /(GAMMA(m)*GAMMA(1+alp)*GAMMA(1+c)));
                    /     m
                    |(E x)  hypergeom([m], [1 + m, 1 - alp], E x)
p := Pi csc(Pi alp) |--------------------------------------------
                    \    GAMMA(c) GAMMA(1 - alp) GAMMA(1 + m)

            c                                      \
       (E x)  hypergeom([c], [1 + c, 1 + alp], E x)|
     - --------------------------------------------|
           GAMMA(m) GAMMA(1 + alp) GAMMA(1 + c)    /

>
> thisp := eval(subs(alp=c-m,p),[c=5,m=1/2,E=25]);
            /     1/2  1/2
            |35 25    x    hypergeom([1/2], [-7/2, 3/2], 25 x)
thisp := Pi |-- ----------------------------------------------
            \64                       Pi

                5                                \
       1562500 x  hypergeom([5], [11/2, 6], 25 x)|
     - ------- ----------------------------------|
         567                   Pi                /

>
> evalf[10](eval(thisp,x=40));
                                              23
                              -0.3141592654 10

> evalf[30](eval(thisp,x=40));
                                      0.

> evalf[50](eval(thisp,x=40));
             0.99999999999999999831055348863866510505212220867081

>
> newthisp:=expand(map(convert,convert(thisp,StandardFunctions),expln)):

> newthisp := simplify(newthisp) assuming x>0;

newthisp := -1/192

           (3/2)         1/2         2                              1/2
    (7000 x      + 1395 x    + 5000 x  + 4350 x + 192 - 192 exp(10 x   ))

             1/2
    exp(-10 x   )

>
> evalf[10](eval(newthisp,x=40));
                                 0.9999999999

> evalf[30](eval(newthisp,x=40));
                       0.999999999999999999999982319406

The following form may make it appear more clear,

> collect(expand(newthisp),exp(x^(1/2))^10);
                           3/2        1/2        2
                      875 x      465 x      625 x    725 x
                    - -------- - -------- - ------ - ----- - 1
                         24         64        24      32
                1 + ------------------------------------------
                                        1/2 10
                                   exp(x   )

A loftier goal might be to get a nice form without having to instantiate at specific values of c, m, and E. Maybe it can be done under assumptions, but I do not yet see how.

acer

Axel Vogt's picture

played with it ...

Comment on 2008-07-17 15:48 from Axel Vogt (Maple Rating 4 885) 
Played with it ... 

In your upload you say 'this is CDF of a distribution' and finally
regarded that ... grrr ...

First set E=1 (as this can be handled through x) and kick off the
leading constant Pi*csc(Pi*a) = Pi/sin(Pi*(c-m)), call that R.

Then diff(R, x) = (cosh(2*X)-sinh(2*X))*polynom(X) / Pi for your
data (where x = X^2), using convert(%,StandardFunctions) as acer
suggests.

This should give a recipe: your pdf is exponential * polynomial.
since exp(-2*X) = (cosh(2*X)-sinh(2*X)) and you should try that
as a starting point, as integrating gives a simple form.

So: may be you can find the pdf distribution, then it is simple.
JacquesC's picture

I agree

Comment on 2008-07-19 09:55 from JacquesC (Maple Rating 4 1648)

In fact, the function might look complicated, but from the right perspective, it is extremely simple:

p := Pi*csc(Pi*alp)*((E*x)^m*hypergeom([m], [1-alp, 1+m], E*x)/(GAMMA(c)*GAMMA(1-alp)*GAMMA(1+m))
-(E*x)^c*hypergeom([c], [1+alp, 1+c], E*x)/(GAMMA(m)*GAMMA(1+alp)*GAMMA(1+c))):
de := PDEtools[dpolyform](y(x)=p, no_Fn);
          3
         d
  de := [--- y(x) =
           3
         dx

                    / 2      \
                    |d       |
        (m + c - 3) |--- y(x)|                         /d      \
                    |  2     |   (m + x - 1 + c - c m) |-- y(x)|
                    \dx      /                         \dx     /
        ---------------------- + -------------------------------]
                  x                             2
                                               x

        &where []

As long as c-m+1/2 is an integer, you will be able to find a nice closed-form for that expression. In other cases, you can reduce things down to the integral of the difference between 2 Bessel functions.

More complex ODE

Comment on 2008-07-19 19:38 from jakubi (Maple Rating 4 689)

I get instead a lengthy 5th order ODE that depends also on alp.

JacquesC's picture

My error

Comment on 2008-07-19 20:56 from JacquesC (Maple Rating 4 1648)

I had alp = m-c defined in my session, as well as having renormalized E*x to x.  Sorry.  With those definitions, you should get the 3rd order ODE.

Comment viewing options

Select your preferred way to display the comments and click "Save settings" to activate your changes.
}