| |
 |
|
|
Science Forum Index » Math - Numerical Analysis Forum » Calculating hermite polynomials
Page 1 of 1
|
| Author |
Message |
| jacob navia |
 Posted: Tue Feb 27, 2007 6:44 am |
|
|
|
Guest
|
Hi people
I am trying to implement the calculation of the hermite polynomials
hermite(n,x)
For x = 1 I get with increasing n this listing. Does this look correct
to you?
Thanks in advance for any help.
[ 0] 1.000000000000000e+000
[ 1] 2.000000000000000e+000
[ 2] 2.000000000000000e+000
[ 3] -4.000000000000000e+000
[ 4] -2.000000000000000e+001
[ 5] -8.000000000000000e+000
[ 6] 1.840000000000000e+002
[ 7] 4.640000000000000e+002
[ 8] -1.648000000000000e+003
[ 9] -1.072000000000000e+004
[ 10] 8.224000000000000e+003
[ 11] 2.308480000000000e+005
[ 12] 2.807680000000000e+005
[ 13] -4.978816000000000e+006
[ 14] -1.725760000000000e+007
[ 15] 1.048916480000000e+008
[ 16] 7.275112960000000e+008
[ 17] -1.901510144000000e+009
[ 18] -2.853840435200000e+010
[ 19] 1.137755648000000e+010
[ 20] 1.107214478336000e+012
[ 21] 1.759326697472000e+012
[ 22] -4.298435469516800e+013
[ 23] -1.633790840791040e+014
[ 24] 1.650522147819520e+015
[ 25] 1.114324033143600e+016
[ 26] -6.023962672810390e+016
[ 27] -6.999277506908820e+017
[ 28] 1.853084341935850e+018
[ 29] 4.290212272256110e+019
[ 30] -2.167464638715710e+019
[ 31] -2.617476656127980e+021
[ 32] -3.891125236252210e+021
[ 33] 1.597362555196860e+023
[ 34] 5.762867766320180e+023
[ 35] -9.709491822074620e+024
[ 36] -5.975905800839050e+025
[ 37] 5.795652951725920e+026
[ 38] 5.581300882966080e+027
[ 39] -3.288436066718480e+028
[ 40] -5.011101902057240e+029
[ 41] 1.628528472963340e+030
[ 42] 4.434809254279600e+031
[ 43] -4.810020664332820e+031
[ 44] -3.910136371967120e+033
[ 45] -3.587454559321360e+033
[ 46] 3.447373643583980e+035
[ 47] 1.019520548174360e+036
[ 48] -3.036627115334070e+037
[ 49] -1.586065149314200e+038
[ 50] 2.658681543164550e+039 |
|
|
| Back to top |
|
| Peter Spellucci |
| Posted: Tue Feb 27, 2007 6:44 am |
|
|
|
Guest
|
In article <45e40b7f$0$5109$ba4acef3@news.orange.fr>,
jacob navia <jacob@jacob.remcomp.fr> writes:
Quote: Hi people
I am trying to implement the calculation of the hermite polynomials
hermite(n,x)
if n==0
p=1;
elseif n==1
p=2*x;
else
p0=1;
p1=2*x;
for i=1:n
p2=2*x*p1-2*i*p0;
p0=p1;
p1=p2;
end
end
Quote:
For x = 1 I get with increasing n this listing. Does this look correct
to you?
Thanks in advance for any help.
o.k.
hth
peter
Quote:
[ 0] 1.000000000000000e+000
[ 1] 2.000000000000000e+000
[ 2] 2.000000000000000e+000
[ 3] -4.000000000000000e+000
[ 4] -2.000000000000000e+001
[ 5] -8.000000000000000e+000
[ 6] 1.840000000000000e+002
[ 7] 4.640000000000000e+002
[ 8] -1.648000000000000e+003
[ 9] -1.072000000000000e+004
[ 10] 8.224000000000000e+003
[ 11] 2.308480000000000e+005
[ 12] 2.807680000000000e+005
[ 13] -4.978816000000000e+006
[ 14] -1.725760000000000e+007
[ 15] 1.048916480000000e+008
[ 16] 7.275112960000000e+008
[ 17] -1.901510144000000e+009
[ 18] -2.853840435200000e+010
[ 19] 1.137755648000000e+010
[ 20] 1.107214478336000e+012
[ 21] 1.759326697472000e+012
[ 22] -4.298435469516800e+013
[ 23] -1.633790840791040e+014
[ 24] 1.650522147819520e+015
[ 25] 1.114324033143600e+016
[ 26] -6.023962672810390e+016
[ 27] -6.999277506908820e+017
[ 28] 1.853084341935850e+018
[ 29] 4.290212272256110e+019
[ 30] -2.167464638715710e+019
[ 31] -2.617476656127980e+021
[ 32] -3.891125236252210e+021
[ 33] 1.597362555196860e+023
[ 34] 5.762867766320180e+023
[ 35] -9.709491822074620e+024
[ 36] -5.975905800839050e+025
[ 37] 5.795652951725920e+026
[ 38] 5.581300882966080e+027
[ 39] -3.288436066718480e+028
[ 40] -5.011101902057240e+029
[ 41] 1.628528472963340e+030
[ 42] 4.434809254279600e+031
[ 43] -4.810020664332820e+031
[ 44] -3.910136371967120e+033
[ 45] -3.587454559321360e+033
[ 46] 3.447373643583980e+035
[ 47] 1.019520548174360e+036
[ 48] -3.036627115334070e+037
[ 49] -1.586065149314200e+038
[ 50] 2.658681543164550e+039 |
|
|
| Back to top |
|
| jacob navia |
| Posted: Tue Feb 27, 2007 9:50 am |
|
|
|
Guest
|
|
| Back to top |
|
| Robert Israel |
| Posted: Tue Feb 27, 2007 1:18 pm |
|
|
|
Guest
|
spellucci@fb04373.mathematik.tu-darmstadt.de (Peter Spellucci) writes:
Quote:
In article <45e40b7f$0$5109$ba4acef3@news.orange.fr>,
jacob navia <jacob@jacob.remcomp.fr> writes:
Hi people
I am trying to implement the calculation of the hermite polynomials
hermite(n,x)
if n==0
p=1;
elseif n==1
p=2*x;
else
p0=1;
p1=2*x;
for i=1:n
p2=2*x*p1-2*i*p0;
p0=p1;
p1=p2;
end
end
For x = 1 I get with increasing n this listing. Does this look correct
to you?
Thanks in advance for any help.
o.k.
hth
peter
[ 0] 1.000000000000000e+000
[ 1] 2.000000000000000e+000
[ 2] 2.000000000000000e+000
[ 3] -4.000000000000000e+000
[ 4] -2.000000000000000e+001
[ 5] -8.000000000000000e+000
[ 6] 1.840000000000000e+002
[ 7] 4.640000000000000e+002
[ 8] -1.648000000000000e+003
[ 9] -1.072000000000000e+004
[ 10] 8.224000000000000e+003
[ 11] 2.308480000000000e+005
[ 12] 2.807680000000000e+005
[ 13] -4.978816000000000e+006
[ 14] -1.725760000000000e+007
[ 15] 1.048916480000000e+008
[ 16] 7.275112960000000e+008
[ 17] -1.901510144000000e+009
[ 18] -2.853840435200000e+010
[ 19] 1.137755648000000e+010
[ 20] 1.107214478336000e+012
[ 21] 1.759326697472000e+012
[ 22] -4.298435469516800e+013
[ 23] -1.633790840791040e+014
[ 24] 1.650522147819520e+015
[ 25] 1.114324033143600e+016
[ 26] -6.023962672810390e+016
[ 27] -6.999277506908820e+017
[ 28] 1.853084341935850e+018
[ 29] 4.290212272256110e+019
[ 30] -2.167464638715710e+019
[ 31] -2.617476656127980e+021
[ 32] -3.891125236252210e+021
[ 33] 1.597362555196860e+023
[ 34] 5.762867766320180e+023
[ 35] -9.709491822074620e+024
[ 36] -5.975905800839050e+025
[ 37] 5.795652951725920e+026
[ 38] 5.581300882966080e+027
[ 39] -3.288436066718480e+028
[ 40] -5.011101902057240e+029
[ 41] 1.628528472963340e+030
[ 42] 4.434809254279600e+031
[ 43] -4.810020664332820e+031
[ 44] -3.910136371967120e+033
[ 45] -3.587454559321360e+033
[ 46] 3.447373643583980e+035
[ 47] 1.019520548174360e+036
[ 48] -3.036627115334070e+037
[ 49] -1.586065149314200e+038
[ 50] 2.658681543164550e+039
For comparison, Maple's values for HermiteH(49,1) and HermiteH(50,1)
to 20 digits are -1.5860651493141992944e+038 and
2.6586815431645459903e+039 respectively, so this seems to be
working quite well.
--
Robert Israel israel@math.MyUniversitysInitials.ca
Department of Mathematics http://www.math.ubc.ca/~israel
University of British Columbia Vancouver, BC, Canada |
|
|
| Back to top |
|
| jacob navia |
| Posted: Tue Feb 27, 2007 1:46 pm |
|
|
|
Guest
|
Robert Israel a écrit :
Quote: spellucci@fb04373.mathematik.tu-darmstadt.de (Peter Spellucci) writes:
In article <45e40b7f$0$5109$ba4acef3@news.orange.fr>,
jacob navia <jacob@jacob.remcomp.fr> writes:
Hi people
I am trying to implement the calculation of the hermite polynomials
hermite(n,x)
if n==0
p=1;
elseif n==1
p=2*x;
else
p0=1;
p1=2*x;
for i=1:n
p2=2*x*p1-2*i*p0;
p0=p1;
p1=p2;
end
end
For x = 1 I get with increasing n this listing. Does this look correct
to you?
Thanks in advance for any help.
o.k.
hth
peter
[snip]
Quote: [ 48] -3.036627115334070e+037
[ 49] -1.586065149314200e+038
[ 50] 2.658681543164550e+039
For comparison, Maple's values for HermiteH(49,1) and HermiteH(50,1)
to 20 digits are -1.5860651493141992944e+038 and
2.6586815431645459903e+039 respectively, so this seems to be
working quite well.
Thanks. I am using a recursive function similar to the one proposed by
Peter Spellucci. The only problem I see is that precision will
eventually be lost when you ask for a <n> of 500, say.
1) Is this actually important? Are there uses of hermite polynomials
where high values of <n> are needed?
2) Are there any other, non-recursive methods?
jacob |
|
|
| Back to top |
|
| Peter Spellucci |
| Posted: Wed Feb 28, 2007 12:53 am |
|
|
|
Guest
|
In article <45e46e82$0$27404$ba4acef3@news.orange.fr>,
jacob navia <jacob@jacob.remcomp.fr> writes:
snip
Quote: Thanks. I am using a recursive function similar to the one proposed by
Peter Spellucci. The only problem I see is that precision will
eventually be lost when you ask for a <n> of 500, say.
1) Is this actually important? Are there uses of hermite polynomials
where high values of <n> are needed?
2) Are there any other, non-recursive methods?
jacob
the error may indeed blow up and this could be important.
I guess that you will not be interested in the values of hermite(n,x) alone
but in sums of the form
sum_n a_n*hermite(n,x)
and if this is the case then you should not compute the polynomial values
separately but use the Clenshaw recursion (quite similar to Horners scheme)
which directly incorporates the recursion for the polynomials into the summation.
Nevertheless such high degree expansions might finally require higher precision.
the presentation of these polynomails in the monomial basis would make things
much worse.
the URL
http://mathworld.wolfram.com/HermitePolynomial.html
has also a "direct" contour integral relation for computing h_n(z)
but I guess this will make things not easier
hth
peter |
|
|
| Back to top |
|
| Kevin G. Rhoads |
| Posted: Wed Feb 28, 2007 10:14 am |
|
|
|
Guest
|
Well, but the contour integral may well be amenable to a steepest-descent
approximation, which may then be implementable either as an analytic
solution in closed form, or some numeric integration. Depends upon the
problem domain and specifics. |
|
|
| Back to top |
|
| Gottfried Helms |
| Posted: Wed Feb 28, 2007 12:35 pm |
|
|
|
Guest
|
Am 27.02.2007 18:46 schrieb jacob navia:
Quote: where high values of <n> are needed?
2) Are there any other, non-recursive methods?
I have a matrix-representation.
Define the simple subdiagonal matrix lGS
lGS =
. . . . . . . .
. . . . . . . .
-1 . . . . . . .
. -3 . . . . . .
. . -6 . . . . .
. . . -10 . . . .
. . . . -15 . . .
. . . . . -21 . .
-----------------------------------------------
then compute the matrix-exponential:
GS = MExp(lGS)
1 . . . . . . .
. 1 . . . . . .
-1 . 1 . . . . .
. -3 . 1 . . . .
3 . -6 . 1 . . .
. 15 . -10 . 1 . .
-15 . 45 . -15 . 1 .
. -105 . 105 . -21 . 1
(This matrix is also known as "modified hermite polynomials"
in mathworld and is directly related to the gaussian normal-function).
Then scale rows and columns by powers of sqrt(2), using the
powerseriesvector V(r) = diag(1,r,r^2,r^3,...) , and r = sqrt(2)
then
HERMITE = V(r) * GS * V(r)
= V(r) * MExp(lGS) * V(r)
= V(r) *( I + lGS/1! + lGS^2/2! .... +lGS^n/n!) * V(r)
-------
In a practical computation the matrixexponential would be
computed by the powerseries expansion, where each term is
again a one-subdiagonal-matrix whose summing is easy to perform.
The V(r) scaling can then be implemented already in the
matrixexpoential-computation by inserting of smooth even powers of r
(that are: powers of 2), since the odd powers vanish because their
cofactors are all zero.
The entries in GS can also be *directly* computed by doublefactorials
and binomial-coefficients, don't have it at hand, but you may visit
my article
http://go.helms-net.de/math/binomial/04_1_gaussmatrix.pdf
where I have given the explicite representation of the entries
of GS.
HTH -
Gottfried Helms |
|
|
| Back to top |
|
| Gottfried Helms |
| Posted: Wed Feb 28, 2007 1:10 pm |
|
|
|
Guest
|
Am 28.02.2007 17:35 schrieb Gottfried Helms:
Quote: Am 27.02.2007 18:46 schrieb jacob navia:
where high values of <n> are needed?
2) Are there any other, non-recursive methods?
(..)
The entries in GS can also be *directly* computed by doublefactorials
and binomial-coefficients, don't have it at hand, but you may visit
my article
http://go.helms-net.de/math/binomial/04_1_gaussmatrix.pdf
where I have given the explicite representation of the entries
of GS.
Well, just looked at my article.
It's as follows:
define for the row- and column-indexes r,c (beginning at zero)
m = (r-c)/2
The unsigned entries of GS and HERMITE are
GS := GS[ r,c] = 0 (if (r-c) is odd)
= (2m-1)!! * binomial(r,c) (if (r-c) is even)
HERMITE:= H[r,c] = GS[r,c]*2^m (if (r-c) is even)
----------------------------------------------------------------
Note, that by using (2n)!/n!/2^n instead of (n-1)!! the
term 2^m cancels out and
the unsigned entries of GS and HERMITE are
GS := GS[ r,c] = 0 (if (r-c) is odd)
= (2m)!/m!/2^m * binomial(r,c) (if (r-c) is even)
HERMITE:= H[r,c] = 0 (if (r-c) is odd)
= (2m)!/m! * binomial(r,c) (if (r-c) is even)
(modulo editing/indexing bugs...)
Gottfried Helms |
|
|
| Back to top |
|
| |
|
Page 1 of 1
All times are GMT - 5 Hours
The time now is Tue Dec 02, 2008 1:01 am
|
|