The following is a specific case of ImportanceSampling.docx in which a = 0 and b = ¥.
Consider the integral
Let
So that
There are a few restrictions on g. First it must be positive definite. Secondly, it must be integrable over an infinite range. The third feature is not strictly necessary, but it would be nice if it were integrable to a finite value in such a way that we can easily find r(t).
I start by looking for a function of x that is an exact derivative,
because then I can integrate it. The
function
Consider
On substituting this into
è
Note that the integral from 0 to
infinity is 1. So that the relation between x and t given by is
Solving for r(t) yields
The r values.doc will use the Laurent transform as its integration scheme. For the forward transform the ψ(r) = exp(-r) is shown in figure 1. For the back transform Y(p) = 1/(1+p2)2 is shown.
The integrand in is finite as tà1 and rè infinity if f(r) approaches zero at least as fast as 1/r2. When f à 0 as exp(-x), GaussLag.doc is a competitive method. The code in exptest.zip is set up to do 5 points
C:\temp\exptest>midpts
NUMERICAL FULL INTEGRAL 0.9991064832699446 – 5 pts
NUMERICAL FULL INTEGRAL 1.0000166663828580 – 50 pts
NUMERICAL FULL INTEGRAL 1.0000001666666860 – 500 pts
NUMERICAL FULL INTEGRAL 1.0000000016666650 – 5000 pts
ANALYTIC FULL INTEGRAL 1.0000000000000000
Figure 1
f(r(t)) versus t for a Laurent transform of
exp(-x)
Figure 2 plot of 1/(1+p(t)2)2
versus t – 20 values. pint.zip
With a slight modification to 20 pts LTpts.zip
NUMERICAL
FULL INTEGRAL 1.0001083036209220
ANALYTIC
FULL INTEGRAL 1.0000000000000000
I R A
1 0.2564102564102564E-01 1.051939513477975
2 0.8108108108108109E-01 1.168736303871439
3 0.1428571428571428 1.306122448979592
4 0.2121212121212122 1.469237832874197
5 0.2903225806451613 1.664932362122789
6 0.3793103448275862 1.902497027348395
7 0.4814814814814815 2.194787379972565
8 0.6000000000000000 2.560000000000000
9
0.7391304347826089 3.024574669187146
10 0.9047619047619050 3.628117913832201
11
1.105263157894737 4.432132963988921
12
1.352941176470589 5.536332179930798
13
1.666666666666667
7.111111111111112
14
2.076923076923078
9.467455621301779
15
2.636363636363638
13.22314049586778
16
3.444444444444445
19.75308641975309
17
4.714285714285716 32.65306122448982
18
7.000000000000000
64.00000000000000
19
12.33333333333334
177.7777777777780
20
39.00000000000014
1600.000000000011
As with Gauss Quadrature
The derivatives at the ends in Figure 1.5 are not equal. This means that Richardson.doc extrapolation can be used.
Figure 3 Point spacing for 2 and 6 point
intervals.
With 21 and 7 points h2-1/21, h1 = 1/7 so
F63 1.0000104981260350
F21
1.0000973605628160
FR
0.9999996403214373
The code for this is in LTint.zip. The code midpts constructs the data set while the cod LTint.for reads the data set, calculates the full sum, the sum with 3*h makes the extrapolation.
IMPLICIT REAL*8(A-H,O-Z)
PARAMETER (N=5000)
DATA PI/3.141592653589793D0/
C
*** MID-POINT
AINT=0
H=1D0/N
open(1,file='ftest.out')
DO I=1,N
T=H*(I-.5D0)
R=T/(1-T)
F=(1+R)**2*FTEST(R)
write(1,'(2g20.6)')t,f
AINT=AINT+F
ENDDO
AINT=AINT*H
ANAL=PI/2
PRINT*,' NUMERICAL FULL INTEGRAL ',AINT
PRINT*,' ANALYTIC FULL INTEGRAL ',ANAL
READ(*,'(A)')
END
To
be specific consider the case
Then
the integral
FUNCTION FTEST(X)
IMPLICIT REAL*8(A-H,O-Z)
FTEST=1/(1+X*X)
RETURN
END
1.570796333461560
ß numerical
1.570796326794900
ß analytical
Figure 4 (1+r(t))*f(r(t)) versus t
Note that the ends are not smooth.
Evaluate both to r and to infinity the integrals
a.
Dwight 856.08
b.