| |
 |
|
|
Science Forum Index » Math - Numerical Analysis Forum » Confluent Hypergeometric/Kummer function calculation
Page 1 of 1
|
| Author |
Message |
| RRogers |
Posted: Mon Mar 19, 2007 9:50 am |
|
|
|
Guest
|
I am solving Lamm's equation and having a problem with Kummer's
functions in a realistic case. Maple does a reasonable (but not
entirely believable) job; the Matlab code fails when run on Octave.
In any case, the series solutions have intrinsic difficulties, and I
am looking atintergal solutions. Unfortunately they also blow up, so
I am trying to repartition the terms via. Mellin xforms.
Is anybody interested enough to help me produce reasonable answers?
My original solutions in Maple look okay, but the person I am talking
with would really like MatLab algorithms. the parameters are:
a= -100.1, -300.1, -1300, -4000 more or less
b=1
z= 5..6
I have acquired a several books; but not all of them.
Ray |
|
|
| Back to top |
|
| user923005 |
Posted: Mon Mar 19, 2007 3:10 pm |
|
|
|
Guest
|
On Mar 19, 7:50 am, "RRogers" <rerog...@plaidheron.com> wrote:
Quote: I am solving Lamm's equation and having a problem with Kummer's
functions in a realistic case. Maple does a reasonable (but not
entirely believable) job; the Matlab code fails when run on Octave.
In any case, the series solutions have intrinsic difficulties, and I
am looking atintergal solutions. Unfortunately they also blow up, so
I am trying to repartition the terms via. Mellin xforms.
Is anybody interested enough to help me produce reasonable answers?
My original solutions in Maple look okay, but the person I am talking
with would really like MatLab algorithms. the parameters are:
a= -100.1, -300.1, -1300, -4000 more or less
b=1
z= 5..6
I have acquired a several books; but not all of them.
The Cephes library by Moshier has 1F1 in it.
/* hyperg.c
*
* Confluent hypergeometric function
*
*
*
* SYNOPSIS:
*
* double a, b, x, y, hyperg();
*
* y = hyperg( a, b, x );
*
*
*
* DESCRIPTION:
*
* Computes the confluent hypergeometric function
*
* 1 2
* a x a(a+1) x
* F ( a,b;x ) = 1 + ---- + --------- + ...
* 1 1 b 1! b(b+1) 2!
*
* Many higher transcendental functions are special cases of
* this power series.
*
* As is evident from the formula, b must not be a negative
* integer or zero unless a is an integer with 0 >= a > b.
*
* The routine attempts both a direct summation of the series
* and an asymptotic expansion. In each case error due to
* roundoff, cancellation, and nonconvergence is estimated.
* The result with smaller estimated error is returned.
*
*
*
* ACCURACY:
*
* Tested at random points (a, b, x), all three variables
* ranging from 0 to 30.
* Relative error:
* arithmetic domain # trials peak rms
* DEC 0,30 2000 1.2e-15 1.3e-16
qtst1:
21800 max = 1.4200E-14 rms = 1.0841E-15 ave = -5.3640E-17
ltstd:
25500 max = 1.2759e-14 rms = 3.7155e-16 ave = 1.5384e-18
* IEEE 0,30 30000 1.8e-14 1.1e-15
*
* Larger errors can be observed when b is near a negative
* integer or zero. Certain combinations of arguments yield
* serious cancellation error in the power series summation
* and also are not in the region of near convergence of the
* asymptotic series. An error message is printed if the
* self-estimated relative error is greater than 1.0e-12.
*
*/ |
|
|
| Back to top |
|
| RRogers |
Posted: Tue Mar 20, 2007 7:51 am |
|
|
|
Guest
|
On Mar 19, 4:10 pm, "user923005" <dcor...@connx.com> wrote:
Quote: On Mar 19, 7:50 am, "RRogers" <rerog...@plaidheron.com> wrote:
I am solving Lamm's equation and having a problem with Kummer's
functions in a realistic case. Maple does a reasonable (but not
entirely believable) job; the Matlab code fails when run on Octave.
In any case, the series solutions have intrinsic difficulties, and I
am looking atintergal solutions. Unfortunately they also blow up, so
I am trying to repartition the terms via. Mellin xforms.
Is anybody interested enough to help me produce reasonable answers?
My original solutions in Maple look okay, but the person I am talking
with would really like MatLab algorithms. the parameters are:
a= -100.1, -300.1, -1300, -4000 more or less
b=1
z= 5..6
I have acquired a several books; but not all of them.
The Cephes library by Moshier has 1F1 in it.
/* hyperg.c
*
* Confluent hypergeometric function
*
*
*
* SYNOPSIS:
*
* double a, b, x, y, hyperg();
*
* y = hyperg( a, b, x );
*
*
*
* DESCRIPTION:
*
* Computes the confluent hypergeometric function
*
* 1 2
* a x a(a+1) x
* F ( a,b;x ) = 1 + ---- + --------- + ...
* 1 1 b 1! b(b+1) 2!
*
* Many higher transcendental functions are special cases of
* this power series.
*
* As is evident from the formula, b must not be a negative
* integer or zero unless a is an integer with 0 >= a > b.
*
* The routine attempts both a direct summation of the series
* and an asymptotic expansion. In each case error due to
* roundoff, cancellation, and nonconvergence is estimated.
* The result with smaller estimated error is returned.
*
*
*
* ACCURACY:
*
* Tested at random points (a, b, x), all three variables
* ranging from 0 to 30.
* Relative error:
* arithmetic domain # trials peak rms
* DEC 0,30 2000 1.2e-15 1.3e-16
qtst1:
21800 max = 1.4200E-14 rms = 1.0841E-15 ave = -5.3640E-17
ltstd:
25500 max = 1.2759e-14 rms = 3.7155e-16 ave = 1.5384e-18
* IEEE 0,30 30000 1.8e-14 1.1e-15
*
* Larger errors can be observed when b is near a negative
* integer or zero. Certain combinations of arguments yield
* serious cancellation error in the power series summation
* and also are not in the region of near convergence of the
* asymptotic series. An error message is printed if the
* self-estimated relative error is greater than 1.0e-12.
*
*/
Thanks! I will look at it. BTW: GNU has a gsl library in C also.
I have sort of given up on summation, and am not familiar enough with
asymptotic expansions to use them without effort. But I am pretty
sure I have resolved the Kummer calculation successfully, and also
resolved the confluent differential equation solutions avoiding most
of the problems and annoyances.
It's not a cookbook approach (which I honestly prefer) but is accurate
and stable.
Ray |
|
|
| Back to top |
|
| |
|
Page 1 of 1
All times are GMT - 5 Hours
The time now is Thu Jan 08, 2009 4:34 am
|
|