|Difficult Integrals: My Original Solutions and Comments|
Message #20 Posted by Valentin Albillo on 11 July 2006, 5:18 p.m.,
in response to message #1 by Valentin Albillo
Thanks for all your interesting inputs, whether actual results or comments, to
my question about the actual efectiveness of HP calcs' built-in methods of integration
when applied to these not-so-easy definite integrals.
I'll give now correct numerical results for all six integrals, and will discuss
what happens when you try to compute them blindly, and what can you do
when this naive approach doesn't deliver. All given code is for the HP-71B, which
uses much the same Romberg's quadrature method as former (HP-34C) and later models,
so your results should be about the same no matter what model you do use.
I1 = | sin(x)*sin(x^2) .dx = 0.491695777984
This is a very difficult integral to compute numerically to any decent accuracy,
because the integrand does not tend to zero as the limit of integration tends to
infinity but oscillates instead with increasing frequency, so at first sight it
might be questionable whether it actually converges or not (which is does).
The naive approach, using 10^6 for infinity, takes a looong time, and blatantly fails:
Insisting on blindly going on by using a plain 4096-point Gaussian quadrature applied to
the interval (0, 10^6), subdivided in 16384 subintervals does succeed in
producing four correct digits, 0.4917, but it takes a long time as well and requires a fast PC, so it seems plain that some math thinking will be needed if we're going to go much farther (or even that far) while using a much slower handheld.
Fortunately, there's a nifty closed form for a generalized family of similar integrals,
this being but one of the simplest cases, namely:
/ A / A
I1 = (Cos(A)* | Cos(x)/Sqrt(x).dx + Sin(A)* | Sin(x)/Sqrt(x).dx )/2
/ 0 / 0
where A=1/4 and both non-elementary integrals are particular cases of Fresnel functions.
Using the above expression, the integral can be easily and accurately computed right from the
keyboard (no programming needed) as follows:
>A=1/4 @ P=1E-12
I2 = | cos(ln(x)/x)/x .dx = 0.323367431678
This one also takes quite long, if using the naive approach, and the initial result is not exactly promising:
Also, this time even a ponderous 8192-point Gaussian formula applied to
4096 subintervals does utterly fail, producing the completely
wrong value 1.7656...
In this case, the easiest way to numerically compute this integral
is by integrating in the complex plane, as discussed previously in
my post of the reference:
"Re: Integration Times "Old" 33s vs "New" 33s
Message #5 Posted by Valentin Albillo on 13 Dec 2005, 4:28 a.m.,
in response to message #1 by John Smitherman"
Self-quoting from that post:
You just need to remember that
e^ix = cos(x) + i*sin(x)
and the original integral becomes:
Integral = Integral(0, 1, RealPart(e^(i*log(x)/x)/x))
= Integral(0, 1, RealPart(x^(i/x-1)))
which can readily be integrated by considering it a line integral in the
complex plane along some integration path. A simple parabolic path will do:
z = t + t*(1-t)*i
which takes us far from the oscillations of the real line.
Using this path, your faithful HP-71B can compute the integral to full
12-digit precision in a few minutes, as follows:
10 DESTROY ALL @ COMPLEX Z @ SFLAG -1 @ DISP INTEGRAL(0,1,1E-9,FNY(IVAR))
20 DEF FNY(T) @ Z=(T,T-T*T) @ FNY=REPT(Z^((0,1)/Z-1)*(1,1-2*T))
which, as stated, is the correct result rounded to the 12 decimal digits shown.
This technique of computing difficult real-line integrals by taking them to the
complex plane is actually very useful in many situations, and is discussed at
length (with an example) in the HP-15C Advanced Functions Handbook.
I3 = | sin(1/x)/x .dx = 0.624713256428
Once again, the naive approach takes quite long, and miserably fails:
Also, the 8192-point Gaussian formula, applied to 1024 subintervals,
fares no better, producing utter garbage: -2.342289...
The way to go is by using 1/x = t, which, together with the well-known elementary
result for the integral between 0 and Infinity of Sin(t)/t .dt, results in the following
/ Inf / Inf / 1 / 1
I3 = | sin(t)/t .dt = | sin(t)/t .dt - | sin(t)/t .dt = PI/2 - | sin(t)/t .dt
/ 1 / 0 / 0 / 0
which can be computed very quickly, to provide 12-digit accuracy:
I4 = | cos(x)*ln(x) .dx = -0.946083070367
For one, the naive approach does work, even if it takes quite long:
Using just one subinterval gets us 10 more-or-less correct digits.
We can get full 12-digit accuracy by subdividing further.
I5 = | ln(1+x)*sin(10*x) .dx = -0.197626807719
The easiest of the pack, the naive approach works too, and very quickly:
and produces nearly 12-digit accuracy.
I6 = | x^(Pi/4)*sin(Pi/4/(2-x) .dx = 1.01123909053
The naive approach also works, but with reduced accuracy and
taking very long:
In search of increased accuracy, we can make it last twice as long by using two subintervals:
>P=PI/4 @ E=1E-12
or even ten times longer:
>S=0 @ FOR X=0 TO 1.9 STEP .1 @ S=S+INTEGRAL(X,X+.1,E,IVAR^P*SIN(P/(2-IVAR))) @ NEXT X @ S
But actually, this is yet another integral which is much easier to compute in the complex plane:
10 DESTROY ALL @ COMPLEX Z @ SFLAG -1 @ DISP INTEGRAL(0,2,1E-10,FNY(IVAR))
20 DEF FNY(T) @ Z=(T,T+T-T*T) @ FNY=REPT((0,-1)*EXP(PI/4*(LOG(Z)+(0,1)/(2-Z)))*(1,2-2*T))
Best regards from V.