On 21 f=A8=A6v, 17:26, Dan Larson <dlar...@[EMAIL PROTECTED]
> wrote:
> On Feb 21, 9:30 am, Spon <christoph.b...@[EMAIL PROTECTED]
> wrote:
>
>
>
> > On Feb 20, 3:44 pm, noahh.schwa...@[EMAIL PROTECTED]
wrote:
>
> > > Hello everyone,
>
> > > I am looking for the ConfluentHypergeometricFunctionof the First
> > > Kind in theIDLMath Library but it does not seem to be implemented!
>
> > > I would like to use afunctionsimilar to the Hypergeometric1F1[a, b,
> > > z] of Mathematica [http://reference.wolfram.com/mathematica/ref/
> > > Hypergeometric1F1.html].
>
> > > I have not found what I was looking for, and so decided to try to
code=
> > > it my self... [sigh...]. Beeing a fresh beginner inIDLthis is a hard
> > > task!
>
> > > Would anybody know how to code an infinite series expansion like the
> > > Hypergeometric1F1?
>
> > > Thank you in advance for your time!
> > > Noah
>
> > Noah,
>
> > here's my attempt. It accepts only scalar inputs for A and B, while Z
> > can be a vector. I've tested it for the examples on the mathematica
> > site and it seems to give correct results, and works correctly for
> > complex input too as far as I can tell. 'Precision' is an input
> > variable to specify how close two successive iterations have to be
> > before thefunctionassumes they are the same and aborts the while
> > loop. Default is 7 (i.e. stop when results differ by 10^-7 or less).
> > If you're finding this programme is running very slow, try decreasing
> > the precision (I was surprised how fast it runs despite the while
> > loop, actually!)
>
> > Ideally the input parameters should all be double precision before you
> > make the call to the funcion, but thefunctionconverts them if
> > they're not.
>
> > If you want all your inputs to be vectors (not just Z), I'm sure it
> > can be done, but it'd be a bit more complicated. :-)
>
> > Take care,
> > Chris
>
> >FUNCTIONHYPERGEOMETRICONEFONE, A, B, Z, $
> > PRECISION =3D Precision, $
> > K =3D K ; K is an output parameter to count No. of WHILE loops
> > performed.
>
> > ; References:
> > ;http://reference.wolfram.com/mathematica/ref/Hypergeometric1F1.html
> > ;http://en.wikipedia.org/wiki/Confluent_hypergeometric_function
>
> > IF N_PARAMS() NE 3 THEN MESSAGE, 'Must input A, B & Z as 3 input
> > parameters.'
> > IF N_ELEMENTS(A) GT 1 THEN MESSAGE, 'Variable A must be a scalar.'
> > IF N_ELEMENTS(B) GT 1 THEN MESSAGE, 'Variable B must be a scalar.'
>
> > A *=3D 1.0D ; Double precision or double complex scalar
> > B *=3D 1.0D ; Double precision or double complex scalar
> > Z *=3D 1.0D ; Double precision or double complex scalar or vector
> > IF N_ELEMENTS(Precision) EQ 0 THEN $
> > Precision =3D 7L ELSE $
> > Precision =3D (LONG(Precision))[0]
> > Cutoff =3D 10D^(-1D * Precision) > (MACHAR()).EPS ; Cutoff can't be
> > smaller than machine accuracy!
>
> > K =3D 0L
> > ThisResult =3D REPLICATE(0D, N_ELEMENTS(Z))
> > WHILE (N_ELEMENTS(LastResult) EQ 0) || (MAX(ABS(LastResult -
> > ThisResult)) GT Cutoff) DO BEGIN
> > LastResult =3D ThisResult
> > AK =3D GAMMA(A + K) / GAMMA(A) ; Define (A)k
> > BK =3D GAMMA(B + K) / GAMMA(B) ; Define (B)k
> > F =3D (AK * Z^K) / (BK * FACTORIAL(K)) ; Evaluatefunction.
> > ThisResult =3D LastResult + F
> > K +=3D 1
> > ENDWHILE ; Until result is good to Precision
>
> > ; Error if not enough while loops to give accurate results.
> > IF K LE 1 THEN MESSAGE, 'Functionfailed. Try greater precision.'
>
> > RETURN, ThisResult
> > END- Hide quoted text -
>
> > - Show quoted text -
>
> Noah,
>
> Here is my implementation, both of the series expansion of
thehypergeometr=
icfunctionand the integral representation. Depending on
> the parameters, I have found that one may be more stable than the
> other. Both of these are based on Arfken and Weber, Mathematical
> Methods for Physicists.
>
> best,
> dan
>
> ; chss
> ; confluenthypergeometricseries solution
> ; calculates the solution to the differential equation:
> ; xy"(x) + (c - x)y'(x) - ay(x) =3D 0
> ; (see Afken and Weber, p. 801-2)
> ;
> ; inputs:
> ; n: number of terms to calculate. Due limits inIDLarchitecture n
> must be between 1 and 170.
> ; a, c: vector of constants - see above. They MUST have the same
> number of elelments
> ; x: input value
> ;
> ; outputs:
> ; y: output vector
> ;
> ; Dan Larson, 2007.10.11
>
> functionchss, n, a, c, x
> y =3D DBLARR((SIZE(a))[1])
> ; first term of series expansion is 1, so we initialize the output
> y[*] =3D 1
>
> FOR i =3D 0, (SIZE(a))[1]-1 DO BEGIN
> ; initialize constant series products
> a_p =3D 1
> c_p =3D 1
> FOR j =3D 0, n DO BEGIN
>
> a_p *=3D (a[i] + j)
> c_p *=3D (c[i] + j)
> d =3D FACTORIAL(j + 1)
>
> y[i] +=3D (a_p/c_p)*(x^(j + 1))/d
>
> ENDFOR
> ENDFOR
> RETURN, y
> END
>
> ; chins
> ; confluenthypergeometricintegral solution
> ; calculates the solution to the differential equation:
> ; xy"(x) + (c - x)y'(x) - ay(x) =3D 0
> ; (see Afken and Weber, p. 801-2)
> ;
> ; inputs:
> ; a, c: vector constants - see above
> ; x: input value
> ;
> ; outputs:
> ; y: output vector
> ;
> ;Dan larson, 2007.10.11
>
> Functionchins, a, c, x
> ; curve point resolution, change this if your results look like
> crap
> b =3D 1000
> ; initialize t vector
> t =3D DINDGEN(b + 1)/b
> ; initialize vector of curve heights
> h =3D DBLARR(b + 1)
> ; initialize output
> y=3Ddblarr(n_elements(c))
>
> g1 =3D GAMMA(c)
> g2 =3D GAMMA(a) * GAMMA(c - a)
>
> FOR i =3D 0, (SIZE(a))[1]-1 DO BEGIN
> FOR j =3D 0, b DO BEGIN
> h[j] =3D exp(x*t[j]) * ((t[j])^(a[i] - 1.0)) * ((1.0 -
> t[j])^(c[i] - a[i] - 1.0))
> ENDFOR
>
> y[i] =3D (g1[i]/g2[i]) * int_tabulated(t, h, /double)
> ENDFOR
>
> RETURN, y
> END
>
> ; chcmp
> ; compares the results between the integral and series form of CHF
> ;
> ; inputs:
> ; a, c: constants
> ; x: input value
> ;
> ; output:
> ; none
> ; (solution is printed to screen)
>
> PRO chcmp, a, c, x, epsilon
> y1 =3D chins(a, c, x)
>
> lastVal =3D 0.0000
>
> FOR n =3D 1, 170 DO BEGIN
> y2 =3D chss( n, a, c, x)
>
> IF ( ABS(y1 - y2)/y1 LT epsilon) THEN BEGIN
> print, "Number of necessary terms for ep =3D ", epsilon, " n
=3D
> ", n
> BREAK
> ENDIF
>
> IF( ABS(lastVal - y2)/y2 EQ 0) THEN BEGIN
> print, "Series converged before matching integral!"
> BREAK
> ENDIF
>
> lastVal =3D y2
> ENDFOR
>
> IF n EQ 169 THEN PRINT, "Equations did not meet!"
>
> RETURN
> END
Thank you all for your answers. It really made my day!
I've tested the 3 methods for the "Mathematica" example (i.e. for
1F1(1,2,[-5..5]))
HYPERGEOMETRICONEFONE, CHSS and CHINS seem to work fine for this
example. HYPERGEOMETRICONEFONE seems a bit faster than the other ones.
Unfortunately I am more interested in evaluating something like 1F1(-.
75,1,40) :
IDL> print,hypergeometriconefone(-0.75D,1D,-40D)
NaN
IDL> print, chins([-0.75D],[1D],-40D)
NaN
IDL> print, chss(169D,[-0.75D],[1D],-40D)
17.478776
The mathematica website gives [http://functions.wolfram.com/
webMathematica/FunctionEvaluation.jsp?name=3DHypergeometric1F1]:
Created by webMathematica =A1=D6 17.5496890473939
Thanks,
Noah
|
