Hello!Currently I am looking into the following integral int_{-infinity}^{infinity} exp{-s^2/2} [ (a*b + s^2)/(a+s^2) ] ^{c} dswhere a and b are positive constants and c is a positive fraction. As long as c is a positive integer I can get to an analytic expression for this integral with a bunch of error functions. However, for fractional values of c I see no analytic expression. My first question is if anybody knows of an analytical expression for this integral. There is a need to evaluate this integral over and over again and computational time starts to play a role. With no analytic expression at hand, my attention has turned to numerical integration. A possible way would be to use an article by Genz and Keister "Fully Symmetric Interpolatory Rules for Multiple Integrals over Infinite Regions with Gaussian Weight" which can be found here There is even a code available in Fortran here However, my knowledge of Fortran is zero and I wonder if somebody has transcribed this algorithm to Matlab / C / VBASo far I work in this one dimensional case with a standard Matlab routine like quadl, but like to see higher precision (that is better than 6 decimals) for testing purposes. The post of Espen Haug in this thread is in my mindgreetsallu

I do not know of an analytic solution, nor any C++ program. However, you should have no trouble converting the Fortran code to VB. There are only minor syntax differences, e.g. DO I = 1, NF DIFFER = ( RESULT(I) - WORK(I) )/CLTOTL WORK(I) = WORK(I) + RNTCLS*DIFFER WORK(NF+I) = ( CLTOTL - RNTCLS )*( WORK(NF+I)/CLTOTL + RNTCLS*DIFFER**2 ) END DOis identical in VB except the first line is:For I = 1 to NFthe last line is:Nextand the ** should be replaced by ^.And you don't have to translate all the code, the real work is done in two fairly short subroutines, the rest is documentation, housekeeping and standard functions.

QuoteOriginally posted by: alluHello!Currently I am looking into the following integral int_{-infinity}^{infinity} exp{-s^2/2} [ (a*b + s^2)/(a+s^2) ] ^{c} dswhere a and b are positive constants and c is a positive fraction. As long as c is a positive integer I can get to an analytic expression for this integral with a bunch of error functions. However, for fractional values of c I see no analytic expression. My first question is if anybody knows of an analytical expression for this integral. There is a need to evaluate this integral over and over again and computational time starts to play a role. With no analytic expression at hand, my attention has turned to numerical integration. A possible way would be to use an article by Genz and Keister "Fully Symmetric Interpolatory Rules for Multiple Integrals over Infinite Regions with Gaussian Weight" which can be found here There is even a code available in Fortran here However, my knowledge of Fortran is zero and I wonder if somebody has transcribed this algorithm to Matlab / C / VBASo far I work in this one dimensional case with a standard Matlab routine like quadl, but like to see higher precision (that is better than 6 decimals) for testing purposes. The post of Espen Haug in this thread is in my mindgreetsalluthere was a French grad student at Berkeley a long time ago that studied the existence of closed form solutions to integrals using Galois theory. trying to remember his name.....

QuoteOriginally posted by: Noniusthere was a French grad student at Berkeley a long time ago that studied the existence of closed form solutions to integrals using Galois theory. trying to remember his name.....Liouville?

the thing here is to avoid taking the whole area of integration at once... the gaussian weight is easily divided into regions where simple asymptotics hold, as is your integrand. You should set up an asymptotic expansion, but to do this you should decide where the variables are going to be, because you have a lot of them. For instance, you can see that the leading term is b^c. When s is large, the tails are easy to estimate, but the cutoff region is going to be a function of a,b,c in order to ensure some uniformity of convergence. Beyond that, there are lots of silly tricks... i.e. use Taylor where applicable, integrate by parts lots of times, etc. Because c is not an integer you may get into some ungly stuff if you try Fourier-type methods, but that's another route since one can find the poles easily and then integrate around a rectangular contour. Since the integrand need only cover the positive reals (by symmetry), maybe you can do some kind of slice-of-pie contour or something...Anyhow, if there is anything interesting to do here, it would come from knowing how to break the infinite region into two pieces... probably somewhere around s = const . c . (log b) ^ (1/2) (i.e. where the gaussian kernel is roughly the size of the integrand).

QuoteOriginally posted by: elanQuoteOriginally posted by: Noniusthere was a French grad student at Berkeley a long time ago that studied the existence of closed form solutions to integrals using Galois theory. trying to remember his name.....Liouville?No, Manuel Brownstein or something like that....

actually, here is the guy....http://www-sop.inria.fr/safir/WHOSWHO/M ... eng.htmlhe was in the office next to mine....

Thanks for the reference Nonius. Indeed his work seems worth reading. Started to look at his course notes which opened an eye without any eureka yet. Interesting enough to see his working place: I was in Sophia Antipolis just a few weeks ago... Anyway, my goal is to have an accurate answer reasonably fast as this integration has to be performed many times for different a, b and c. Thinking along the way KR suggested: take a bounded interval from 0 to a critical value of s. The remaining can be shown to be small in my case, not dependent on a,b,c. The remaining integral is bounded but not solvedGuess the easiest way out is to go to the Fortran code Genz is providing. Yes, the Fortran structure is the same as VBA though initialization is worrysome. Thinking about transcribing the whole thing-- allu

QuoteOriginally posted by: alluThanks for the reference Nonius. Indeed his work seems worth reading. Started to look at his course notes which opened an eye without any eureka yet. Interesting enough to see his working place: I was in Sophia Antipolis just a few weeks ago... Anyway, my goal is to have an accurate answer reasonably fast as this integration has to be performed many times for different a, b and c. Thinking along the way KR suggested: take a bounded interval from 0 to a critical value of s. The remaining can be shown to be small in my case, not dependent on a,b,c. The remaining integral is bounded but not solvedGuess the easiest way out is to go to the Fortran code Genz is providing. Yes, the Fortran structure is the same as VBA though initialization is worrysome. Thinking about transcribing the whole thing-- alluI couldn't get to the notes, what is the gist of his research.....it sounded pretty bizarre the first time I heard about it...but hey, he got a Springer book published.

My thoughts: with Maple i do not have problems to get numerical values, so itshould not be a problem with Matlab.To do it by 'manual coding' a change of variables x/sqrt(2)=xi gives youInt(exp(-xi^2)*((a*b+2*xi^2)/(a+2*xi^2))^c*2^(1/2),xi = -infinity .. infinity)For that you find a recipe in Abramowitz&Stegun (25.4.46) by a Gauss-Hermiteif you restrict to a finite range -bound ... +bound. The bounds here dependon your a and b. As a finite sum the method that is easy to code, the neccessary zeros of the Hermites you find in Abramowitz (up to degree 20)and you find an error estimation as well (somewhat ugly in this case).---Edited after reading it again:Gauss-Hermite is not the best idea, applying Simpson's rule is easy and works.Writing int(f(x), x=0..infinity) ~ int(f(x), x=A..B), A=0, B=A+n*h withstepsize h and n even(!) steps as approximating sum according to that rule reads as1/3*h*(f(A)+f(A+n*h)+4*Sum(f(A+(2*i-1)*h),i = 1 .. 1/2*n)+2*Sum(f(A+2*h*i),i = 1 .. 1/2*n-1))Checking on an example: for a=1,b=2,c=3/2 take the upper bound B=8 and n=64. Then2.695974507757100 resp 2.695974507845684 resp 2.695974507845685 are the values forSimpson resp int(..., x=0..8) resp int(..., x=0..infinity) using 16 digits exactness.

I am still interested in knowning whether a closed form solution exists for the integral. Please advise.

I don't see any good contours, even with integral c. The problem is the rapid decay of the exponential term - the function has to blow up if you try to slide the contour past a pole to a different latitude - basically I don't see any good trade-in candidates for the original contour (well, like I said below, I just used the +ve reals by symmetry). This is one of those moments when I wish I owned a copy of Erdelyi et al.... but I don't. Besides that there is always brute-force Taylor, which is not a bad idea and will give you good convergence as a function of the 'a' parameter. Write (...)^c as 1 + (b-1).(1/(1+s^2/a))... you may have to fiddle with the 'b' variable in order to get something nice when you apply taylor to f(x) = (b' + x)^c. Then use geom series for the part that decays like 1/s^2. You will ultimately arrive at integrals of powers of s against the gaussian kernel, which have a closed form. The powers will be increasing, but b/c of the gaussian kernel the individual integrals are not increasing quickly in size. However, each one will give you a higher power of 'a' in the denominator...

Thanks people for all the comments. The solution is getting closer!For integral c the problem seems ok and Maple shows you some answers without too much hesitation. For your info here is c=1: 1/2*a/b*2^(1/2)*(-1/2*Pi*sqrt(b)*sqrt(2)*exp(1/2*b)*erf(1/2*sqrt(2)*sqrt(b))+1/2*Pi*sqrt(2)*sqrt(b)*exp(1/2*b))+1/2*sqrt(b)*Pi^(1/2)*2^(1/2)/b~+Pi*exp(1/2*b)*erf(1/2*sqrt(2)*sqrt(b))-Pi*exp(1/2*b)and c=21/2*a^2/b^2*2^(1/2)*(1/2*sqrt(Pi)*b+1/12*Pi*sqrt(b)*(-3+3*b)*sqrt(2)*exp(1/2*b)*erf(1/2*sqrt(2)*sqrt(b))+1/4*Pi*sqrt(2)*sqrt(b)*(1-b)*exp(1/2*b))+2*a/b^2*2^(1/2)*(-1/4*sqrt(Pi)*b^2-1/8*Pi*b^(3/2)*sqrt(2)*(1+b)*exp(1/2*b)*erf(1/2*sqrt(2)*sqrt(b))+1/8*Pi*b^(3/2)*sqrt(2)*(1+b)*exp(1/2*b))+2/b^2*2^(1/2)*(-3/16*Pi*b^(5/2)*sqrt(2)*(1+1/3*b)*exp(1/2*b)+1/4*sqrt(Pi)*b^2+1/8*sqrt(Pi)*b^3+1/16*Pi*b^(5/2)*sqrt(2)*(3+b)*exp(1/2*b)*erf(1/2*sqrt(2)*sqrt(b)))etc etc the expressions just get larger but remain explicit in error function. Didn't check the methodology for integral c. So far my experience showed Maple is almost always correct; just have to verify! Bronstein might work. His tutorial was 35 pages. Looks readable. AVt: Genz and Keister used Gauss-Hermite rules... Before I used quad (based on adaptive Simpson) and quadl (based on adaptive Lobatto) in Matlab for an extended version of the shown integral. At that moment the routine seemed unstable. Will have to check if it stayed within my error estimation. The question is ofcourse: was the first guess not as bad as I thought -- allu

I do not see why the antiderivative of the integrand should be of a nice closed form.Calling it g it seems at least to satisfy the diff equation g''(x) = -x * g'(x) *( 1+ 2*c*a*(b-1)/((a+x^2)*(a*b+x^2)) ) . The 1st part stands for the cdf normal andthe 2nd can be expressed as a reasonable taylor series (alternating and only in odddegrees). Solving for summands results in powers + hypergeometric fcts so i gave up.May be integral transforms for the DE are better. It seems to be 2nd order Liouvilleand those are seldom nice :-)Allu: i am not in numerics, but then you have several solution and Genz seems onlyneccessary for the multivariate case.