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
Hi all:
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.
/ Inf
I1 = | sin(x)*sin(x^2) .dx = 0.491695777984
/ 0
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:
INTEGRAL(0,1E6,1E-12,SIN(IVAR)*SIN(IVAR*IVAR))
-.356719258
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
>(COS(A)*INTEGRAL(0,A,P,COS(IVAR)/SQR(IVAR))+SIN(A)*INTEGRAL(0,A,P,SIN(IVAR)/SQR(IVAR)))/2
.491695777984
/ 1
I2 = | cos(ln(x)/x)/x .dx = 0.323367431678
/ 0
This one also takes quite long, if using the naive approach, and the initial result is not exactly promising:
>INTEGRAL(0,1,1E-12,COS(LN(IVAR)/IVAR)/IVAR)
.169114784125
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))
>RUN
.323367431678
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.
/ 1
I3 = | sin(1/x)/x .dx = 0.624713256428
/ 0
Once again, the naive approach takes quite long, and miserably fails:
>INTEGRAL(0,1,1E-12,SIN(1/IVAR)/IVAR)
-1.6116828454
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
expression:
/ 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:
>PI/2-INTEGRAL(0,1,1E-12,SIN(IVAR)/IVAR)
.624713256433
/ 1
I4 = | cos(x)*ln(x) .dx = -0.946083070367
/ 0
For one, the naive approach does work, even if it takes quite long:
>INTEGRAL(0,1,1E-12,COS(IVAR)*LN(IVAR))
-.946083070191
Using just one subinterval gets us 10 more-or-less correct digits.
We can get full 12-digit accuracy by subdividing further.
/ 2*Pi
I5 = | ln(1+x)*sin(10*x) .dx = -0.197626807719
/ 0
The easiest of the pack, the naive approach works too, and very quickly:
>INTEGRAL(0,2*PI,1E-12,LN(1+IVAR)*SIN(10*IVAR))
-.197626807716
and produces nearly 12-digit accuracy.
/ 2
I6 = | x^(Pi/4)*sin(Pi/4/(2-x) .dx = 1.01123909053
/ 0
The naive approach also works, but with reduced accuracy and
taking very long:
>INTEGRAL(0,2,1E-12,IVAR^(PI/4)*SIN(PI/4/(2-IVAR)))
1.0112556592
In search of increased accuracy, we can make it last twice as long by using two subintervals:
>P=PI/4 @ E=1E-12
>INTEGRAL(0,1,E,IVAR^P*SIN(P/(2-IVAR)))+INTEGRAL(1,2,E,IVAR^P*SIN(P/(2-IVAR)))
1.01119507067
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
1.01122722759
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))
>RUN
1.01123909053
Best regards from V.
|