The Museum of HP Calculators

HP Forum Archive 20

[ Return to Index | Top of Index ]

Question about Numerical Integration
Message #1 Posted by Namir on 10 Nov 2011, 5:53 p.m.

The Integrate function in several HP calculators calculates the integral between two values, say A and B. Are there any tricks to use the same method when integrating between A and infinity, minus infinity and A, and minus infinity and plus infinity?

I know there are special types of Gaussian quadrature methods that can handle the above cases. My question is directed at working with the algorithms for finite integrals.

Namir

      
Re: Question about Numerical Integration
Message #2 Posted by Valentin Albillo on 10 Nov 2011, 6:10 p.m.,
in response to message #1 by Namir

Quote:
The Integrate function in several HP calculators calculates the integral between two values, say A and B. Are there any tricks to use the same method when integrating between A and infinity, minus infinity and A, and minus infinity and plus infinity?

I know there are special types of Gaussian quadrature methods that can handle the above cases. My question is directed at working with the algorithms for finite integrals.

Namir

The usual trick is to first perform a simple change of variable which will reduce any interval, including infinities at one or both extremes, to any finite interval you care for such as [0,1] or [-1,1].

Best regards from V.

            
Re: Question about Numerical Integration
Message #3 Posted by Lyuka on 10 Nov 2011, 6:41 p.m.,
in response to message #2 by Valentin Albillo

Many RF engineers would like the conversion below,

s = (z - 1) / (z + 1)

that is used to plot various impedance in a Smith chart.

Lyuka

                  
Re: Question about Numerical Integration
Message #4 Posted by Namir on 10 Nov 2011, 11:31 p.m.,
in response to message #3 by Lyuka

looks like:

s = (x^2 - 1) / (x^2 + 1)

is a better transformation, since it is valid for all real values of x.

Namir

      
Re: Question about Numerical Integration
Message #5 Posted by peacecalc on 10 Nov 2011, 11:58 p.m.,
in response to message #1 by Namir

Hallo namir,

when you need only numerical results, you can also use instead of infinities great/small numbers like +/- 10^6.

That's sometimes tricky, shown in the advanced handbook for the 15c.

sincerely peacecalc

Edited: 10 Nov 2011, 11:59 p.m.

            
Re: Question about Numerical Integration
Message #6 Posted by Namir on 11 Nov 2011, 3:15 a.m.,
in response to message #5 by peacecalc

That's the kind of tricks I was looking for.

Putting it in pseudo-code form in the case of integrating from A to infinity:

Given f(x), A, SmallValue, RelTolerance, and DiffTolerance
InfVal=10^6
IntegralVal1 = Integral of f(x) from A to InfVal
Do
  InfVal = 10 * InfVal
  IntegralVal2 = Integral of f(x) from A to InfVal
  if |IntegralVa2| < SmallValue then
    RelativeErr = 0
    DiffErr = IntegralVa2 - IntegralVal1
  else
    RelativeErr = (IntegralVal2 - IntegralVal1) / IntegralVal1
    DiffErr = 0
  end
  IntegralVal1 = IntegralVal2
Until |RelativeErr| < RelTolerance OR |DiffErr| < DiffTolerance
Integral = IntegralVal1

Here is a perhaps more efficient version that uses integration by parts:

Given f(x), A, SmallValue, RelTolerance, and DiffTolerance
B = 10^6 // or any other different high value (10^3, 10^4, 10^5, 10^6, 10^7, etc
IntegralVal0 = Integral of f(x) from A to B // should calculate most of the final answer
InfVal= 10 * B
IntegralVal1 = Integral of f(x) from B to InfVal
Do
  InfVal = 10 * InfVal
  IntegralVal2 = Integral of f(x) from A to InfVal
  if |IntegralVa2| < SmallValue then
    RelativeErr = 0
    DiffErr = IntegralVa2 - IntegralVal1
  else
    RelativeErr = (IntegralVal2 - IntegralVal1) / IntegralVal1
    DiffErr = 0
  end
  IntegralVal1 = IntegralVal2
Until |RelativeErr| < RelTolerance OR |DiffErr| < DiffTolerance
Integral = IntegralVal0 + IntegralVal1

Edited: 11 Nov 2011, 5:47 a.m.

                  
Re: Question about Numerical Integration
Message #7 Posted by Mike (Stgt) on 11 Nov 2011, 7:04 a.m.,
in response to message #6 by Namir

IIRC, there is a discussion worth to read in the PPC-ROM manual.

Hope this hepls

Ciao.....Mike

                  
Re: Question about Numerical Integration
Message #8 Posted by Gjermund Skailand on 11 Nov 2011, 7:11 a.m.,
in response to message #6 by Namir

Tanh-Sinh transformation or "Double exponential method" may also be efficient, especially when the function is oscillating. There is a fast implementation in C (hpgcc2) for HP50G in my package INT1D at www.hpcalc.org.

      
Re: Question about Numerical Integration
Message #9 Posted by Dieter on 11 Nov 2011, 7:53 a.m.,
in response to message #1 by Namir

As already pointed out, there are basically two approaches for handling improper integrals:

  • Use a large (finite) value where the integrand becomes negligible compared to the final result.

    Example: for an 8-digit result of the left tail Normal CDF a lower limit of -6 usually is sufficient as the PDF here is mererly 6E-9 and the whole integral from -infinity to -6 is less than 1E-9.

  • Apply a variable transformation to make infinity "less infinite". ;-) This is the preferred method suggested in most of the other posts in this thread.

    Example: Take a look at the 15C Advanced Functions Handbook, p. 54 ff. There, this topic is discussed with several examples on how to use the 15C's Integrate function for improper integrals. On p. 60-64 you will find a very nice and illustrative application of these ideas, evaluating the Normal CDF as well as the error function, with both plus or minus infinity as the limits, giving exact results even far out in the tails. This is accomplished by using the transformation u = exp(t²) as soon as the integration limit exceeds a certain value (here: 1,6).

Dieter

[ Return to Index | Top of Index ]

Go back to the main exhibit hall