Integrator Manual


 Overview
 
 
 
 
XROM  Function  Desciption
 16,00
 16,01
 16,02
 16,03
 16,04
 16,05
 16,06
 16,07
 16,08
 16,09
 16,10
 16,11
 16,12
 16,13
 16,14
 16,15
 16,16
 16,17
 16,18
 16,19
 16,20
 16,21
 16,22
 16,23
 16,24
 16,25
 16,26
 16,27
 16,28
 16,29
 16,30
 16,31
 16,32
 16,33
 16,34
 16,35
 16,36
 16,37
 16,38
 16,39
 16,40
 16,41
 16,42
 16,43
 16,44
 16,45
 16,46
 16,47
 16,48
 16,49
 16,50
 16,51
 16,52
 16,53
 16,54
 16,55
 16,56
 16,57
 16,58
 16,59
 16,60
 16,61
 16,62
 16,63
 -INTEGRATOR
 "FILON"
 "GCHB1"
 "GCHB2"
 "GH20"
 "GH20D"
 "GH20T"
 "GH30"
 "GH30D"
 "GH30T"
 "GK7"
 "GL3"
 "GL3D"
 "GL3M"
 "GL3T"
 "GL10"
 "GL10D"
 "GL10T"
 "GLA10"
 "GLA10D"
 "GLA10T"
 "GLA15"
 "GLA15D"
 "GLA15T"
  GW
  GX
 -DISCRETE
 "3DLS"
 "DIGD"
 "DIGT"
 "ITGCUB"
 "LAGR"
 "LAGRI"
 "NCSI"
 "SIMP"
 "TRAP"
 -ELLIPTIC
 "CEI"
 "ELI"
  ELIPF
 "LEI1"
 "LEI2"
 "LEI3"
  RF
  RFZ
 "RG"
  RJ
  RJZ
 -ROMBERG
 "CURVE"
 "IG"
 "IG3"
 "IGD"
 "IGT"
 "ROM"
 "SREVL"
-HP & PPC
 "DIFEQ"
 "IGAB"
 "IGAI"
 "ITG"
 "PPCIG"
 "PPCSV"
 "RGDATA"
 Section Header
 Filon Integration
 Gauss-Chebyshev Integrals 1st kind
 Gauss-Chebyshev Integrals 2nd kind
 Gauss-Hermite 20-point formula
 Gauss-Hermite 20-point formula - Double Integrals
 Gauss-Hermite 20-point formula - Triple Integrals
 Gauss-Hermite 30-point formula
 Gauss-Hermite 30-point formula - Double Integrals
 Gauss-Hermite 30-point formula - Triple Integrals
 Gauss-Kronrod 7-point ( & 15-point ) formulae
 Gauss-Legendre 3-point formula
 Gauss-Legendre 3-point formula - Double Integrals
 Gauss-Legendre 3-point formula - Multiple Integrals
 Gauss-Legendre 3-point formula - Triple Integrals
 Gauss-Legendre 10-point formula
 Gauss-Legendre 10-point formula - Double Integrals
 Gauss-Legendre 10-point formula - Triple Integrals
 Gauss-Laguerre 10-point formula
 Gauss-Laguerre 10-point formula - Double Integrals
 Gauss-Laguerre 10-point formula - Triple Integrals
 Gauss-Laguerre 15-point formula
 Gauss-Laguerre 15-point formula - Double Integrals
 Gauss-Laguerre 15-point formula - Triple Integrals
 Weights for Gauss-Legendre 10-point formula
 Abscissas for Gauss-Legendre 10-point formula
 Section Header  also loads data constants
 Tridiagonal Linear Systems
 Double Integrals ( equally-spaced arguments )
 Triple Integrals ( equally-spaced arguments )
 Integrals ( connecting cubic segments )
 Lagrange Interpolation Formula
 Integral of the Lagrange Polynomial
 Natural Cubic Spline Integration
 Simpson rule ( unequally-spaced arguments too )
 Trapezoidal rule ( unequally-spaced arguments too )
 Section Header  Assign it to a key...
 Complete Elliptic Integrals ( 1st & 2nd kinds )
 Incomplete Elliptic Integrals ( 1st & 2nd kinds )
 Elliptic Integral of the 1st kind
 Legendre Elliptic Integral of the 1st kind
 Legendre Elliptic Integral of the 2nd kind
 Legendre Elliptic Integral of the 3rd kind
 Carlson Elliptic Integral ( 1st kind )
 Carlson Elliptic Integral ( 1st kind ) - complex arg
 Symmetric Carlson Elliptic Integral of the 2nd kind
 Carlson Elliptic Integral ( 3rd kind )
 Carlson Elliptic Integral ( 3rd kind ) - complex arg
 Section Header  Also stores Z-Y-X into R01-R02-R03
 Arc length of a curve
 Romberg Integration - midpoint formula
 Romberg Integration - midpoint formula - step tripl.
 Midpoint formula - Double Integrals
 Midpoint formula - Triple Integrals
 A subroutine called by many programs in this section
 Area of a surface of revolution
 Section Header
 Differential Equations
 Integral from a to b
 Integral from a to infinity
 Integral
 PPC ROM Integrator
 PPC Solver
 Stores Data Points for discrete cases.

 
 

Constants for Gauss-Laguerre 10-point formula
 
 
 
  •  R07 = 3.084411158 E-01       •  R08 = 1.377934705 E-01
  •  R09 = 4.011199292 E-01       •  R10 = 7.294545495 E-01
  •  R11 = 2.180682876 E-01       •  R12 = 1.808342902
  •  R13 = 6.208745610 E-02       •  R14 = 3.401433698
  •  R15 = 9.501516975 E-03       •  R16 = 5.552496140
  •  R17 = 7.530083886 E-04       •  R18 = 8.330152747
  •  R19 = 2.825923350 E-05       •  R20 = 11.84378584
  •  R21 = 4.249313985 E-07       •  R22 = 16.27925783
  •  R23 = 1.839564824 E-09       •  R24 = 21.99658581
  •  R25 = 9.911827220 E-13       •  R26 = 29.92069701

 

Constants for Gauss-Laguerre 15-point formula
 
 
 
  •  R07 = 2.182348859 E-01     •  R08 = 0.09330781202
  •  R09 = 3.422101779 E-01     •  R10 = 0.4926917403
  •  R11 = 2.630275779 E-01     •  R12 = 1.215595412
  •  R13 = 1.264258181 E-01     •  R14 = 2.269949526
  •  R15 = 4.020686492 E-02     •  R16 = 3.667622722
  •  R17 = 8.563877804 E-03     •  R18 = 5.425336627
  •  R19 = 1.212436147 E-03     •  R20 = 7.565916227
  •  R21 = 1.116743923 E-04     •  R22 = 10.12022857
  •  R23 = 6.459926762 E-06     •  R24 = 13.13028248
  •  R25 = 2.226316907 E-07     •  R26 = 16.65440771
  •  R27 = 4.227430385 E-09     •  R28 = 20.77647890
  •  R29 = 3.921897267 E-11     •  R30 = 25.62389423
  •  R31 = 1.456515264 E-13     •  R32 = 31.40751917
  •  R33 = 1.483027051 E-16     •  R34 = 38.53068331
  •  R35 = 1.600594906 E-20     •  R36 = 48.02608557

 

Constants for Gauss-Hermite 20-point formula
 
 
 
  •  R08 = 4.622436696 E-01      •  R09 = 0.2453407083
  •  R10 = 2.866755054 E-01      •  R11 = 0.7374737285
  •  R12 = 1.090172060 E-01      •  R13 = 1.234076215
  •  R14 = 2.481052089 E-02      •  R15 = 1.738537712
  •  R16 = 3.243773342 E-03      •  R17 = 2.254974002
  •  R18 = 2.283386360 E-04      •  R19 = 2.788806058
  •  R20 = 7.802556479 E-06      •  R21 = 3.347854567
  •  R22 = 1.086069371 E-07      •  R23 = 3.944764040
  •  R24 = 4.399340992 E-10     •  R25 = 4.603682450
  •  R26 = 2.229393646 E-13      •  R27 = 5.387480890

 

Constants for Gauss-Hermite 30-point formula
 
 
 
  •  R08 = 3.863948895 E-01     •  R09 = 0.2011285765
  •  R10 = 2.801309308 E-01     •  R11 = 0.6039210586
  •  R12 = 1.467358475 E-01     •  R13 = 1.008338271
  •  R14 = 5.514417687 E-02     •  R15 = 1.415527800
  •  R16 = 1.470382970 E-02     •  R17 = 1.826741144
  •  R18 = 2.737922473 E-03     •  R19 = 2.243391468
  •  R20 = 3.483101243 E-04     •  R21 = 2.667132125
  •  R22 = 2.938725229 E-05     •  R23 = 3.099970530
  •  R24 = 1.579094887 E-06     •  R25 = 3.544443873
  •  R26 = 5.108522451 E-08     •  R27 = 4.003908604
  •  R28 = 9.178580424 E-10     •  R29 = 4.483055357
  •  R30 = 8.106186297 E-12     •  R31 = 4.988918969
  •  R32 = 2.878607081 E-14     •  R33 = 5.533147152
  •  R34 = 2.810333603 E-17     •  R35 = 6.138279220
  •  R36 = 2.908254700 E-21     •  R37 = 6.863345294

 

Filon Integration
 

-This method computes integrals of the form:    §ab  f(x).cos(k.x).dx   and   §ab  f(x).sin(k.x).dx

    §ab  f(x).cos(k.x).dx  ~  h.[ A(k.h) (  f(x2n).sin(k.x2n)  -  f(x0).sin(k.x0) ) + B(k.h) C2n + C(k.h) C2n-1 ]                 a = x0 ;  b = x2n ;  h = (b-a)/(2n)
    §ab  f(x).sin(k.x).dx   ~  h.[ A(k.h) ( -f(x2n).cos(k.x2n) + f(x0).cos(k.x0) ) + B(k.h) S2n + C(k.h) S2n-1 ]

where   C2n = f(x0).cos(k.x0) + f(x2).cos(k.x2) + ......... + f(x2n).cos(k.x2n)  -  (1/2) ( f(x2n).cos(k.x2n) + f(x0).cos(k.x0) )
             S2n = f(x0).sin(k.x0) + f(x2).sin(k.x2) + ......... + f(x2n).sin(k.x2n)  -  (1/2) ( f(x2n).sin(k.x2n) + f(x0).sin(k.x0) )

            C2n-1 = f(x1).cos(k.x1) + f(x3).cos(k.x3) + ......... + f(x2n-1).cos(k.x2n-1)
            S2n-1 = f(x1).sin(k.x1) + f(x3).sin(k.x3) + ......... + f(x2n-1).sin(k.x2n-1)

 and     A(µ) = 1/µ + (sin(2µ))/(2µ2) - 2(sin2µ)/µ3
           B(µ) =  2(1+cos2µ)/µ2 - 2(sin(2µ))/µ3
           C(µ) =  4(sinµ)/µ3 - 4(cosµ)/µ2

-If  k = 0 , Simpson's rule is used
 

Data Registers:         R00 = Function name

                                    R01 = a
                                    R02 = b         R06 = k             R07 to R12: temp
                                    R03 = n

                          >>>  At the end,    R04 =   §ab  f(x).cos(k.x).dx  = X-register  ,  R05 =   §ab  f(x).sin(k.x).dx   = Y-register
Flags:  /
Subroutine:   A program which computes  f(x)  assuming x is in X-register upon entry.
 
 
         STACK            INPUTS         OUTPUTS
          Alpha           F name                /
              T                a                /
              Z                b                /
              Y                n   §ab  f(x).sin(k.x).dx
              X                k   §ab  f(x).cos(k.x).dx

 
Example:   Evaluate   I = §16  Ln(x).cos(10x).dx    and   J = §16  Ln(x).sin(10x).dx
 
 
 01  LBL "T"
 02  LN
 03  RTN

 
 alpha  "T"  alpha

  1    ENTER^
  6    ENTER^
  8    ENTER^
 10   XEQ "FILON"     >>>>    I =  -0.047890755
                                    X<>Y   J =   0.175512930

with n = 16  ,  16  ENTER^
                      10      R/S     >>>>    I =  -0.047429223
                                          X<>Y   J =   0.174731804

with n = 32  ,  32  ENTER^
                      10      R/S     >>>>    I =  -0.047453034
                                          X<>Y   J =   0.174714501

with n = 64  ,  64  ENTER^
                      10      R/S     >>>>    I =  -0.047454443     ( exact results are  I = -0.047454534
                                          X<>Y   J =   0.174713854                           and  J =  0.174713817   to 9 decimals )
 

Gauss-Chebyshev Integrals 1st kind
 

-This program calculates the singular integral:    §ab  ((x-a)(x-b))-1/2  f(x).dx  ~  (pi/n) [ f(x1) + f(x2) + ..... + f(xn) ]

    where   xi = (a+b)/2 + ((b-a)/2).cos((2i-1)pi/(2n))
 

Data Registers:          R00 = Function name

                                    R01 = a
                                    R02 = b
                                    R03 = n          R04 = Integral         R05 thru R06: temp
Flags:  /
Subroutine:   A program which computes  f(x)  assuming x is in X-register upon entry.
 
 
      STACK       INPUTS     OUTPUTS
        alpha        f name            /
           Z            a            /
           Y            b            /
           X            n       Integral

 
Example:     Compute  §13  ex(-x2+4.x-3)-1/2 dx  =  §13  ex((3-x)(x-1))-1/2 dx
 
 
 01  LBL "T"
 02  E^X
 03  RTN

 
 alpha   "T"  alpha           and, if  n = 2

  1    ENTER^
  3    ENTER^
  2    XEQ "GCHB1"  >>>>  29.262

    n = 4       4          R/S             >>>>  29.389695
    n = 8       8          R/S             >>>>  29.38969917
   n = 16     16         R/S             >>>>  29.38969918

Note:

-This program works in all angular mode,  but "GCHB1" and the subroutine must be executed in the same angular mode.
 

Gauss-Chebyshev Integrals 2nd kind
 

-This program calculates the singular integral:    §ab  ((x-a)(x-b)) 1/2  f(x).dx  ~  [ w1 f(x1) + w2 f(x2) + ..... + wn f(xn) ]

    where   xi = (a+b)/2 + ((b-a)/2).cos((i.pi/n)   and  wi = (pi/(n+1)) Sin2 ((i.pi/n)
 

Data Registers:           R00 = Function name

                                    R01 = a
                                    R02 = b
                                    R03 = n          R04 = Integral         R05 thru R07: temp
Flags:  /
Subroutine:   A program which computes  f(x)  assuming x is in X-register upon entry.
 
 
      STACK       INPUTS     OUTPUTS
        alpha        f name            /
           Z            a            /
           Y            b            /
           X            n       Integral

 
Example:     Compute  §13  ex(-x2+4.x-3) 1/2 dx  =  §13  ex((3-x)(x-1)) 1/2 dx
 
 
 01  LBL "T"
 02  E^X
 03  RTN

 
 alpha   "T"  alpha           and, if  n = 2

  1    ENTER^
  3    ENTER^
  2    XEQ "GCHB2"  >>>>  13.088

    n = 4       4        R/S             >>>>  13.11926565
    n = 8       8        R/S             >>>>  13.11926681
   n = 16     16       R/S             >>>>  13.11926681

Note:

-This program works in all angular mode,  but "GCHB2" and the subroutine must be executed in the same angular mode.
 

Gauss-Hermite 20-point formula
 

-Gauss-Hermite m-point formulas are:    §-¥+¥ exp(-x2) f(x).dx  ~  w1.f(x1) + w2.f(x2) + ....... + wm.f(xm)
  where wi and xi  are choosen so that the formula produces perfect accuracy when  f(x) is a polynomial of degree < 2m
 

Data Registers:             R00 = function name

                                        R01 = Integral                                      R02-R03: temp   R08 to R27 = Constants

Flags: F00 to F07
Subroutine:  A program that takes x in X-register and returns f(x) in X-register
 
 
       STACK        INPUT      OUTPUT
         Alpha       F name           /
            X            /            I

   where   I = §-¥+¥ exp(-x2) f(x).dx

Example:   Evaluate   I = §-¥+¥  exp(-x2)  Ln ( 1 + x + x2 )  dx
 
 
  LBL "T"
  ENTER^
  X^2
  +
  LN1+X
  RTN

 
 Alpha  "T"  Alpha    XEQ "GH20"   >>>>   I ~  0.451490093
 

Gauss-Hermite 20-point formula, Double Integrals
 

-The Gauss-Hermite 20-point & 30-point formulae may also be used to compute   I = §-¥+¥ §-¥+¥ exp(-x2-y2) f(x,y) dx dy
 

Data Registers:             R00 = function name

                                        R01 = Integral                                      R02 to R05: temp   R08 to R27 = Constants
 

Flags: F00 to F07
Subroutine:  A program that takes x in X-register & y in Y-register and returns f(x,y) in X-register
 
 
       STACK        INPUT      OUTPUT
         Alpha       F name           /
            X            /            I

   where   I = §-¥+¥ §-¥+¥ exp(-x2-y2) f(x,y) dx dy

Example:   Evaluate    I = §-¥+¥ §-¥+¥exp(-x2-y2)  Cos ( x + y ) dx dy
 
 
  LBL "T"
  +
  COS
  RTN

 
  XEQ "RAD"    Alpha "T" Alpha    XEQ "GH20D"   >>>>   I ~  1.905472266
 

Gauss-Hermite 20-point formula, Triple Integrals
 

-We try to evaluate    I = §-¥+¥ §-¥+¥ §-¥+¥ exp(-x2-y2-z2) f(x,y,z) dx dy dz    with the 20-point formula
 

Data Registers:             R00 = function name

                                        R01 = Integral                                      R02 to R07: temp   R08 to R27 = Constants

Flags: F00 to F07
Subroutine:  A program that takes x in X-register , y in Y-register & z in Z-register and returns f(x,y,z) in X-register
 
 
       STACK        INPUT      OUTPUT
         Alpha       F name           /
            X            /            I

   where   I = §-¥+¥ §-¥+¥ §-¥+¥ exp(-x2-y2-z2) f(x,y,z) dx dy dz

Example:   Evaluate   I = I = §-¥+¥ §-¥+¥ §-¥+¥ exp(-x2-y2-z2)   Cos ( x + y + z ) dx dy dz
 
 
  LBL "T"
  +
  +
  COS
  RTN

 
  XEQ "RAD"   Alpha "T" Alpha    XEQ "GH20T"   >>>>   I ~  2.630291901

-The exact result is  I = ( PI / sqrt(e) )3/2 = 2.630291900   ( rounded to 9 decimals )

Notes:

-I've used V41 to get this result.
-With a real HP-41, the execution time is probably of the order of 2 hours.
 

Gauss-Hermite 30-point formula
 

-We compute   I  =  §-¥+¥ exp(-x2) f(x).dx
 

Data Registers:             R00 = function name

                                        R01 = Integral                                      R02-R03: temp   R08 thru R37 = Constants

Flags: F00 to F07
Subroutine:  A program that takes x in X-register and returns f(x) in X-register
 
 
       STACK        INPUT      OUTPUT
         Alpha       F name           /
            X            /            I

   where   I = §-¥+¥ exp(-x2) f(x).dx

Example:   Evaluate   I = §-¥+¥  exp(-x2)  Ln ( 1 + x + x2 )  dx
 
 
  LBL "T"
  ENTER^
  X^2
  +
  LN1+X
  RTN

 
  Alpha  "T"  Alpha    XEQ "GH30"   >>>>   I ~  0.451471191

Note:

-With the change of argument  x = u / 1.5 ,  "GH20"  gives  0.451469515  and  "GH30"  returns  0.451469593  which is correct to 9D.
 

Gauss-Hermite 30-point formula, Double Integrals
 

  I = §-¥+¥ §-¥+¥ exp(-x2-y2) f(x,y) dx dy
 

Data Registers:             R00 = function name

                                        R01 = Integral                                      R02 to R05: temp   R08 thru R37 = Constants

Flags: F00 to F07
Subroutine:  A program that takes x in X-register & y in Y-register and returns f(x,y) in X-register
 
 
       STACK        INPUT      OUTPUT
         Alpha       F name           /
            X            /            I

   where   I = §-¥+¥ §-¥+¥ exp(-x2-y2) f(x,y) dx dy

Example:   Evaluate again    I = §-¥+¥ §-¥+¥ exp(-x2-y2)  Cos ( x + y ) dx dy
 
 
  LBL "T"
  +
  COS
  RTN

 
  XEQ "RAD"    Alpha  "T"  Alpha    XEQ "GH30D"   >>>>   I ~  1.905472265
 

-"GH30D"  gives the exact result =  PI / sqrt(e) = 1.905472265  rounded to 9 decimals.
 

Gauss-Hermite 30-point formula, Triple Integrals
 

    I = §-¥+¥ §-¥+¥ §-¥+¥ exp(-x2-y2-z2) f(x,y,z) dx dy dz
 

Data Registers:             R00 = function name

                                        R01 = Integral                                      R02 to R07: temp   R08 thru R37 = Constants
 

Flags: F00 to F07
Subroutine:  A program that takes x in X-register , y in Y-register & z in Z-register and returns f(x,y,z) in X-register
 
 
       STACK        INPUT      OUTPUT
         Alpha       F name           /
            X            /            I

   where   I = §-¥+¥ §-¥+¥ §-¥+¥ exp(-x2-y2-z2) f(x,y,z) dx dy dz

Example:   Evaluate   I = I = §-¥+¥ §-¥+¥ §-¥+¥ exp(-x2-y2-z2)   Cos ( x + y + z ) dx dy dz
 
 
  LBL "T"
  +
  +
  COS
  RTN

 
  XEQ "RAD"   Alpha   "T"  Alpha    XEQ "GH30T"   >>>>   I ~  2.630291900

-The exact result is  I = ( PI / sqrt(e) )3/2 = 2.630291900   ( rounded to 9 decimals )

Notes:

-I've again used V41 to get this result.
-With a real HP-41, the execution time is probably more than 6 hours.
 

Gauss 7-point formula, embedded in a Kronrod 15-point formula
 

-When we use a Gauss-Legendre m-point formula to evaluate I = §ab f(x) dx , there is no error estimate.
-To get an error estimate, we can subdivide the interval into 2 parts and take the difference between the 2 results as an error-estimate.

-Gauss-Kronrod integration adds new abscissas to those already employed by a Gauss-Legendre m-point formula to build a more accurate formula
-Thus, less function-evaluations are required to get an error estimate.

-The following programs use a 7-point Gauss formula embedded in a 15-point Kronrod formula.
-So, 15 function evaluations are computed instead of 21.
 

-The interval [a;b] may be divided into  n sub-intervals  [ a+k.(b-a)/n ; a+(k+1).(b-a)/n ]    k = 0 , 1 , ...... , n-1
-  n must be stored in R03
 

Data Registers:            R00 = function name

                                       R01 = a                                R04 = Gauss7-Integral                                  R06 to R25: temp
                                       R02 = b                                R05 = Kronrod15-Integral
                                       R03 = n
Flags: /
Subroutine:  A program that returns f(x) in X-register, assuming x is in X-register upon entry.
 
 
        STACK        INPUTS      OUTPUTS
        Alpha         F name           /
           Z             a           /
           Y             b           I15
           X             n            I7

  with       I7 = Gauss7-Integral = R04
  and       I15 = Kronrod15-Integral = R05

Example:    The same one: compute  §04 Sin ( x2 ) dx
 
 
  LBL "T"
  X^2
  SIN
  RTN
  END

 
 Alpha  "T"  Alpha            XEQ "RAD"

   0    ENTER^
   4    ENTER^                                                            and if you choose  n = 2
   2    XEQ "GK7"  >>>>   I7 = 0.747111478
                              RDN   I15 = 0.747133845

-With  n = 4

 4      R/S          >>>>   I7 = 0.747133835
                         RDN   I15 = 0.747133845
 

Gauss-Legendre 3-point formula
 

-Gauss-Legendre formulas are defined by  §-11 f(x).dx  ~  w1.f(x1) + w2.f(x2) + ....... + wm.f(xm)
  where wi and xi  are choosen so that the formula produces perfect accuracy when  f(x) is a polynomial of degree < 2m

-A linear change of variable maps an interval [a,b] to the standard [-1,+1]
-The interval [a,b] may also be divided into n subintervals for a better accuracy:
-Theoretically, the approximations tend to the exact result as n tends to infinity.

-"GL3" evaluates  §ab  f(x) dx  with the 3-point formula.
 

Data Registers:         R00 = Function name

                                   R01 = a
                                   R02 = b
                                   R03 = n = number of subintervals over which the formula is applied.          R04 =  §ab f(x).dx        R05 thru R08: temp
Flags:  /
Subroutine:   A program which computes  f(x)  assuming x is in X-register upon entry.
 
 
        STACK        INPUTS      OUTPUTS
        Alpha         F name             /
           Z             a             /
           Y             b             /
           X             n             I

    where   I = §ab  f(x) dx

Example:     Evaluate  §13  e-x^2 dx

-Load this short routine:
 
 
 01  LBL "T" 
 02  X^2
 03  CHS
 04  E^X
 05  RTN

 
  Alpha  "T"  Alpha

  1    ENTER^
  3    ENTER^                                                                     and with  n = 2
  2    XEQ "GL3"  >>>>  0.139390854

       n = 4       4      R/S          >>>>  0.139383255
       n = 8       8      R/S          >>>>  0.139383216

-If a or b is infinite, make a change of variable ( like x = Tan u ) to reduce [ a ; b ] to a finite interval.
 

Gauss-Legendre 3-point formula, Double Integrals
 

   I = §ab §u(x)v(x)  f(x;y) dx dy
 

-The result is in X-register and in R08
 

REGISTERS:    R00 thru R18
FLAG:              F03
 

         R00 = the name of the subroutine which calculates f(x;y)
         R01 = a
         R02 = b
         R03 = n = the number of subintervals into which the intervals (a;b) and (u(x);v(x)) are to be divided.
    •   R04 = the name of the subroutine which calculates u(x).
    •   R05 = the name of the subroutine which calculates v(x).

 >>>  3 subroutines must be keyed into program memory:

      -The first one must take x and y from the X-register and Y-register respectively and calculate f(x;y).
      -The second takes x from the X-register and calculates u(x).
      -The third takes x from the X-register and calculates v(x).
 
 
        STACK        INPUTS      OUTPUTS
        Alpha         F name             /
           Z             a             /
           Y             b             /
           X             n             I

    where   I =  §ab §u(x)v(x)  f(x;y) dx dy

Example:   Evaluate  I = §12 §xx^2 (1+x4 y4)1/2 dx dy.

We must load the 3 needed subroutines into program memory ( any global labels, maximum of 6 characters ), for instance:
 
 
  01  LBL "FF"
  02  *
  03  X^2
  04  X^2
  05  SIGN
  06  LAST X
  07  +
  08  SQRT
  09  RTN
  10  LBL "U"
  11  RTN
  12  LBL "V"
  13  X^2
  14  RTN

 
     "U"  ASTO 04
     "V"  ASTO 05

  "FF"  Alpha                   and with n = 1

   1    ENTER^
   2    ENTER^
   1    XEQ "GL3D"  the result is  15.45937082    in the X-register and in R08

             n = 2     2     R/S    >>>>   15.46673275
             n = 4     4     R/S    >>>>   15.46686031
             n = 8     8     R/S    >>>>   15.46686245  ... all the digits are correct!
 

Gauss-Legendre 3-point formula, Triple Integrals
 

   I  =  §ab §u(x)v(x) §w(x;y)t(x;y)   f (x;y;z) dx dy dz
 

-The result is in X-register and in R08
 

REGISTERS:    R00 thru R18
FLAG:              F03
 

         R00 = the name of the subroutine which calculates f(x;y;z)
         R01 = a
         R02 = b
         R03 = n = the number of subintervals into which the intervals of integration are to be divided.
    •   R04 = the name of the subroutine which calculates u(x).
    •   R05 = -------------------------------------------  v(x).
    •   R06 = ------------------------------------------- w(x;y)
    •   R07 = ------------------------------------------- t(x;y)

>>>  5 subroutines must be keyed into program memory:

      -The first one must take x, y and z from the X-register,the Y-register and the Z-register respectively and calculate f(x;y;z).
      -The second takes x from the X-register and calculates u(x).
      -The third takes x from the X-register and calculates v(x).
      -Another one takes x and y from the X and Y registers respectively and calculates w(x;y)
      -The last one -------------------------------------------------------------------- t(x;y)
 
 
        STACK        INPUTS      OUTPUTS
        Alpha         F name             /
           Z             a             /
           Y             b             /
           X             n             I

    where   I = §ab §u(x)v(x) §w(x;y)t(x;y)   f (x;y;z) dx dy dz

Example:      Evaluate   I = §12 §xx^2 §x+yxy  xyz  ( x2 + y2 + z2 ) -1/2 dx dy dz

  The 5 required subroutines are, for instance:
  ( with global labels of 6 characters maximum )
 
 
  01  LBL "FF"
  02  ENTER^
  03  X^2
  04  R^
  05  ST* Z
  06  X^2
  07  +
  08  R^
  09  ST* Z
  10  X^2
  11  +
  12  SQRT
  13  /
  14  RTN
  15  LBL "U"
  16  RTN
  17  LBL "V"
  18  X^2
  19  RTN
  20  LBL "W"
  21  +
  22  RTN
  23  LBL "T"
  24  *
  25  RTN

 
     "U"    ASTO 04
     "V"    ASTO 05
    "W"    ASTO 06
     "T"    ASTO 07

     "FF"  Alpha            and with n = 1

       1    ENTER^
       2    ENTER^
       1    XEQ "GL3T"  gives    I1 =  0.765014888

         n = 2      2           R/S     >>>>        I2 =  0.770640690
         n = 4      4           R/S     >>>>        I4 =  0.770731245
         n = 8      8           R/S     >>>>        I8 =  0.770732669

Note:

-Registers R23 , R24 , .... may be employed to compute
 the functions if they are too much complicated to fit in the stack.
 

Gauss-Legendre 3-point formula, Multiple Integrals
 

-"GL3M" calculates:    I = §ab §u1(x1)v1(x1)  §u2(x1;x2)v2(x1;x2) .........  §u(n-1)(x1;...;x(n-1))v(n-1)(x1;...;x(n-1))   f(x1;x2;....;xn)  dx1dx2....dxn   ( n < 7 )

   i-e integrals , double integrals , ...... , up to sextuple integrals, provided the subroutines do not contain any XEQ.
-This limitation is only due to the fact that the return stack can accomodate up to 6 pending return addresses.
-The 3 point Gauss-Legendre formula is applied to n sub-intervals in the direction of all axis:  §-11 f(x).dx = ( 5.f(-(3/5)1/2) + 8.f(0) + 5.f((3/5)1/2) )/9
-Linear changes of arguments transform the intervals into the standard interval  [ 1 ; 1 ]
 
 

Data Registers:       R00 = m =  number of  §    ;   R01 = x1 ;  R02 = x2 ; .......... ;  R06 = x6
                          R07 = 7 ; R08 = 4 ; R09 = 2 ; R10 = 1.6 ; R11 = 3.6 ; R12 = 0.151/2 ; R13 = 0.61/2   ;    R14 thru R17 = control numbers

                     R18 = n = number of sub-intervals

                     R19 = f name                                  R25 = u1 name         R32 = u2 name         .....................    R53 = u5 name
                     R20 = a                                          R26 = v1 name         R33 = v2 name        .....................     R54 = v5 name
                     R21 = b                                          R27 = u1(x1)            R34 = u2(x1;x2)       .....................     R55 = u5(x1;x2;...;x5)
                     R22 = (b-a)/n                                 R28 = v1(x1)            R35 = v2(x1;x2)       .....................      R56 = v5(x1;x2;...;x5)
                     R23 = index                                   R29 = (v1-u1)/n       R36 = (v2-u2)/n        .....................     R57 = (v5-u5)/n
                     R24 = partial sum,                          R30 = index             R37 = index             .....................      R58 = index
               and, finally, the integral                          R31 = partial sum     R38 = partial sum    ......................      R59 = partial sum
 

Flags: /
Subroutines:  The functions  f ;  u1 ; v1 ; u2 ; v2 ; ....... are to be computed assuming  x1 is in R01 ; x2 is in R02 ; ........
 
 
       STACK        INPUT      OUTPUT
            X            /            I

 
Example:    Evaluate   I  =  §13  §xx^2  §x+yx.y  §zx+z   ln(x2+y/z+t) dx.dy.dz.dt

-First, we load the 7 required subroutines ( I've used one-character global labels to speed up execution, namely "T" , "U" , "V" , "W" , "X" , "Y" , "Z" ):
 
 
  LBL "T"             f(x,y,z,t) = ln(x2+y/z+t)
  RCL 01
  X^2
  RCL 02
  RCL 03
  /
  +
  RCL 04
  +
  LN
  RTN
  LBL "U"           u1(x) = x
  RCL 01
  RTN
  LBL "V"           v1(x) = x2
  RCL 01
  X^2
  RTN
  LBL "W"          u2(x,y) = x+y
  RCL 01
  RCL 02
  +
  RTN
  LBL "X"          v2(x,y) = x.y
  RCL 01
  RCL 02
  *
  RTN
  LBL "Y"          u3(x,y,z) = z
  RCL 03
  RTN
  LBL "Z"          v3(x,y,z) = x+z
  RCL 01
  RCL 03
  +
  RTN

 
-We execute SIZE 046 ( or greater ) and

      XEQ "GL3M"   >>>>    " A=?"
         1       R/S       >>>>     " B=?"
         3       R/S       >>>>     " FNAME?"
         T       R/S      >>>>      " U1NAME?"
         U       R/S     >>>>      " V1NAME?"
         V       R/S     >>>>      " U2NAME?"
         W      R/S     >>>>      " V2NAME?"
         X       R/S     >>>>      " U3NAME?"
         Y       R/S     >>>>      " V3NAME?"
         Z        R/S     >>>>      " U4NAME?"     press  R/S   without any entry when all function names have been keyed in
                   R/S     >>>>      " N=?"                ( number of sub-intervals )  let's try     n = 1

         1        R/S     >>>>      I = 160.4523151     in X-register and in register R24

-To recalculate I with another n-value, key in this value and press XEQ 10 or simply press R/S , the message "N=?" will be displayed. Enter it and R/S ...

  for example,

    with  n = 2     2  XEQ 10  yields   160.6314953        ( or  R/S  >>>  "N=?"   2   R/S   >>>>   160.6314953 )
    and   n = 4     4  XEQ 10  yields   160.6342731

Notes:

-If n is multiplied by 2, execution time is approximately multiplied by 16 for quadruple integrals , by 64 for sextuple integrals.
-We can use Lagrange interpolation formula to improve our results by extrapolation to the limit.
-In this example, we find  I = 160.634317
 

Gauss-Legendre 10-point formula
 

-Gauss-Legendre formulas are defined by  §-11 f(x).dx  ~  w1.f(x1) + w2.f(x2) + ....... + wm.f(xm)
  where wi and xi  are choosen so that the formula produces perfect accuracy when  f(x) is a polynomial of degree < 2m

-A linear change of variable maps an interval [a,b] to the standard [-1,+1]
-The interval [a,b] may also be divided into n subintervals for a better accuracy:
-Theoretically, the approximations tend to the exact result as n tends to infinity.

-"GL10" evaluates  §ab  f(x) dx  with the 10-point formula.
 
 

Data Registers:             R00 = function name

                                        R01 = a                              R04 = Integral
                                        R02 = b
                                        R03 = n                              R05 thru R09: temp
Flags: F01 & F02
Subroutine:  A program that takes x in X-register and returns f(x) in X-register
 
 
        STACK        INPUTS      OUTPUTS
        Alpha         F name             /
           Z             a             /
           Y             b             /
           X             n             I

    where   I = §ab  f(x) dx

Example:   Evaluate  §04  Sin x2 dx
 
 
   LBL "T" 
   X^2
   SIN
   RTN

 
   XEQ "RAD"

  Alpha  "T"  Alpha               and if you choose  n = 1

   0   ENTER^
   4   ENTER^
   1   XEQ "GL10"  >>>>   I = 0.748650151

-With n = 2           2          R/S          >>>>   I = 0.747133892

-With n = 3           3          R/S          >>>>   I = 0.747133845       ( correct to 9D )
 

Gauss-Legendre 10-point formula, Double Integrals
 

-"GL10D" uses the Gauss-Legendre 10-point formula  to evaluate:   §ab §u(x)v(x)  f(x;y) dx dy
 

Data Registers:             R00 = f name             ( Registers R04-R05 are to be initialized before executing "GL10D" )

                                        R01 = a                      •  R04 = u name                            R07 thru R18: temp
                                        R02 = b                      •  R05 = v name
                                        R03 = n                          R06 = Integral
Flags: F01 & F02
Subroutines:

-A program that takes x in X-register & y in Y-register and returns f(x,y) in X-register
-A program that takes x in X-register and returns  u(x)  in X-register
-A program that takes x in X-register and returns  v(x)  in X-register
 
 
        STACK        INPUTS      OUTPUTS
        Alpha         F name             /
           Z             a             /
           Y             b             /
           X             n             I

  where  I = §ab §u(x)v(x)  f(x;y) dx dy

Example:   Evaluate  I = §12 §xx^2 (1+x4 y4)1/2 dx dy.
 
 
 01  LBL "T" 
 02  *
 03  X^2
 04  X^2
 05  SIGN
 06  LAST X
 07  +
 08  SQRT
 09  RTN
 10  LBL "U"
 11  RTN
 12  LBL "V"
 13  X^2
 14  RTN

 
  "U"  ASTO 04
  "V"  ASTO 05

   "T"  Alpha              and with n = 1

    1   ENTER^
    2   ENTER^
    1   XEQ "GL10D"  >>>>   15.46686246    in X-register & R06

-With    n = 2        2          R/S           >>>>   15.46686247     in X-register & R06
 

Gauss-Legendre 10-point formula, Triple Integrals
 

-"GL10T" uses the Gauss-Legendre 10-point formula  to evaluate:   §ab §u(x)v(x)  §w(x,y)t(x,y)  f(x,y,z) dx dy dz
 

Data Registers:             R00 = f name            ( Registers R04 thru R07 are to be initialized before executing "GL10T" )

                                        R01 = a                      •  R04 = u name                      •  R06 = w name
                                        R02 = b                      •  R05 = v name                      •  R07 = t name
                                        R03 = n                          R08 = Integral                         R09 to R27: temp
Flags: F01 & F02
Subroutines:

-A program that takes x , y , z  in registers X , Y , Z and returns f(x,y,z) in X-register
-A program that takes x in X-register and returns  u(x)  in X-register
-A program that takes x in X-register and returns  v(x)  in X-register
-A program that takes x , y  in registers X & Y and returns  w(x,y)  in X-register
-A program that takes x , y  in registers X & Y and returns  t(x,y)  in X-register
 
 
        STACK        INPUTS      OUTPUTS
        Alpha         F name             /
           Z             a             /
           Y             b             /
           X             n             I

  where  I = §ab §u(x)v(x)  §w(x,y)t(x,y)  f(x,y,z) dx dy dz

Example:   Evaluate  I = §12 §xx^2 §x+yx.y ( x2+y2+z2 )1/2 dx dy dz
 
 
 01  LBL "S"
 02  X^2
 03  X<>Y
 04  X^2
 05  +
 06  X<>Y
 07  X^2
 08  +
 09  SQRT
 10  RTN
 11  LBL "U"
 12  RTN
 13  LBL "V"
 14  X^2
 15  RTN
 16  LBL "W"
 17  +
 18  RTN
 19  LBL "T"
 20  *
 21  RTN

 
  "U"  ASTO 04
  "V"  ASTO 05

  "W" ASTO 06
  "T"  ASTO 07

   "S"  Alpha              and with n = 1

    1   ENTER^
    2   ENTER^
    1   XEQ "GL10T"  >>>>   0.813302277    in X-register & R08

    n = 2         2        R/S           >>>>   0.813302283     in X-register & R08
    n = 4         4        R/S           >>>>   0.813302278     in X-register & R08

Note:

-Due to roundoff-errors, it's difficult to decide which result is best ! may be the first one ?
 

Gauss-Laguerre 10-point formula
 

-Gauss-Laguerre m-point formulas are of the type:    §0+¥ exp(-x) f(x).dx  ~  w1.f(x1) + w2.f(x2) + ....... + wm.f(xm)
  where wi and xi  are choosen so that the formula produces perfect accuracy when  f(x) is a polynomial of degree < 2m

-"GLA10" employs the 10-point formula.
 

Data Registers:             R00 = function name

                                        R01 = Integral                R02: temp     R07 thru R26 = Constants

Flags: F00 to F07
Subroutine:  A program that takes x in X-register and returns f(x) in X-register
 
 
       STACK         INPUT      OUTPUT
        Alpha         F name            /
            X             /            I

   where   I = §0+¥  exp(-x) f(x) dx

Example:   Evaluate   I = §0+¥  exp(-x)  Ln ( 1 + x )  dx
 
 
  LBL "T"
  LN1+X
  RTN

 
 Alpha  "T"  Alpha    XEQ "GLA10"   >>>>   I ~  0.596354677

Note:

-If you have to calculate §a+¥  exp(-x) f(x) dx , use the change of variable  x = a + u
 

Gauss-Laguerre 10-point formula, Double Integrals
 

   I = §0+¥ §0+¥  exp(-x-y) f(x,y) dx dy
 

Data Registers:             R00 = function name

                                        R01 = Integral                                      R02 to R04: temp     R07 thru R26 = Constants

Flags: F00 to F07
Subroutine:  A program that takes x in X-register & y in Y-register and returns f(x,y) in X-register
 
 
       STACK         INPUT      OUTPUT
        Alpha         F name            /
            X             /            I

   where   I = §0+¥ §0+¥  exp(-x-y) f(x,y) dx dy

Example:   Evaluate   I = §0+¥ §0+¥  exp(-x-y)  Sin ( x + y ) dx dy
 
 
  LBL "T"
  +
  SIN
  RTN

 
  XEQ "RAD"    Alpha  "T"  Alpha    XEQ "GLA10D"   >>>>   I ~  0.500000715
 

Gauss-Laguerre 10-point formula, Triple Integrals
 

   I = §0+¥ §0+¥  §0+¥  exp(-x-y-z) f(x,y,z) dx dy dz
 

Data Registers:             R00 = function name

                                        R01 = Integral                              R02 to R06: temp  R07 thru R26 = Constants

Flag: F01
Subroutine:  A program that takes x in X-register , y in Y-register , z in Z-register and returns f(x,y,z) in X-register
 
 
       STACK         INPUT      OUTPUT
        Alpha         F name            /
            X             /            I

   where   I = §0+¥ §0+¥  §0+¥  exp(-x-y-z) f(x,y,z) dx dy dz

Example:   Evaluate   I = §0+¥§0+¥  §0+¥  exp(-x-y-z)  Sin ( x + y + z ) dx dy dz
 
 
  LBL "T"
  +
  +
  SIN
  RTN

 
  XEQ "RAD"    Alpha  "T"  Alpha    XEQ "GLA10T"   >>>>   I ~  0.250000765
 

Gauss-Laguerre 15-point formula
 

  I  =  §0+¥ exp(-x) f(x).dx
 

Data Registers:             R00 = function name

                                        R01 = Integral                                      R02: temp   R07 thru R36 = Constants

Flags:  F00 to F07
Subroutine:  A program that takes x in X-register and returns  f(x) in X-register
 
 
       STACK         INPUT      OUTPUT
        Alpha         F name            /
            X             /            I

   where   I = §0+¥  exp(-x) f(x) dx

Example:   Evaluate      I = §0+¥  exp(-x)  Ln ( 1 + x )  dx
 
 
  LBL "T"
  LN1+X
  RTN

 
 Alpha  "T"  Alpha    XEQ "GLA15"   >>>>   I ~  0.596347721
 

Notes:

-Is it the best evaluation of this integral that we can get from these formulae ?
-Let's try the change of variable x = u/2.1  the subroutine becomes:
 
 
  LBL "T"
  2.1
  /
  ENTER^
  LN1+X
  X<>Y
  1.1
  *
  E^X
  *
  2.1
  /
  RTN

 
  XEQ "GLA15"  >>>>  I ~ 0.5963473625   whereas  "GLA10"  gives  0.596347379

-The exact result - rounded to 10D - is  0.5963473623
-This trick will work sometimes ... but not always !
 

Gauss-Laguerre 15-point formula, Double Integrals
 

   I = §0+¥ §0+¥  exp(-x-y) f(x,y) dx dy
 

Data Registers:             R00 = function name

                                        R01 = Integral                                      R02 to R04: temp   R07 thru R36 = Constants

Flags:  F00 to F07
Subroutine:  A program that takes x in X-register & y in Y-register and returns f(x,y) in X-register
 
 
       STACK         INPUT      OUTPUT
        Alpha         F name            /
            X             /            I

   where   I = §0+¥ §0+¥  exp(-x-y) f(x,y) dx dy

Example:   Evaluate again    I = §0+¥ §0+¥  exp(-x-y)  Sin ( x + y ) dx dy
 
 
  LBL "T"
  +
  SIN
  RTN

 
  XEQ "RAD"   Alpha  "T"  Alpha    XEQ "GLA15D"   >>>>   I ~  0.5000000001
 

Gauss-Laguerre 15-point formula, Triple Integrals
 

   I = §0+¥ §0+¥  §0+¥  exp(-x-y-z) f(x,y,z) dx dy dz
 

Data Registers:             R00 = function name

                                         R01 = Integral                                      R02 to R06: temp   R07 thru R36 = Constants

Flags:  F00 to F07
Subroutine:  A program that takes x in X-register , y in Y-register , z in Z-register and returns f(x,y,z) in X-register
 
 
       STACK         INPUT      OUTPUT
        Alpha         F name            /
            X             /            I

   where   I = §0+¥ §0+¥  §0+¥  exp(-x-y-z) f(x,y,z) dx dy dz

Example:   Evaluate   I = §0+¥§0+¥  §0+¥  exp(-x-y-z)  Sin ( x + y + z ) dx dy dz
 
 
  LBL "T"
  +
  +
  SIN
  RTN

 
  XEQ "RAD"   Alpha "T" Alpha    XEQ "GLA15T"   >>>>   I ~  0.2500000000

Note:

-This is the exact result: no roundoff-error at all !
 

Weights for Gauss-Legendre 10-point formula
 

 -The M-code function "GW" takes an integer between 1 and 5 and replaces it by
 
 
  X = 1  ----->   X = 0.2955242247
  X = 2  ----->   X = 0.2692667193
  X = 3  ----->   X = 0.2190863625
  X = 4  ----->   X = 0.1494513492
  X = 5  ----->   X = 0.06667134431

 

Abscissas for Gauss-Legendre 10-point formula
 

-The M-code function "GX" takes an integer between 1 and 5 and replaces it by
 
 
  X = 1  ----->   X = 0.1488743390
  X = 2  ----->   X = 0.4333953941
  X = 3  ----->   X = 0.6794095683
  X = 4  ----->   X = 0.8650633667
  X = 5  ----->   X = 0.9739065285

 

Tridiagonal Linear Systems
 

-"3DLS" is called as a subroutine by "NCSI"
-The coefficients are to be stored in column order, and the system is solved without pivoting

REGISTERS:    R00 & Rbb thru Ree ,  ee =  4n + bbb - 3
FLAGS:            /
 
 
       STACK         INPUTS       OUTPUTS
            X     (bbb.eee)system    (bbb.eee)solution

   where  (bbb.eee)system  is the control number of the system    eee =  4n + bbb - 3  , bb > 00 ,  n = number of unknowns ,  n > 1
              (bbb.eee)solution  --------------------------  solution

Example:

   2.x + 5.y                                     = 2                                                  2   5                           2                 R01   R03                                            R17
   3.x + 7.y + 4.z                            = 4                                                  3   7   4                      4                 R02   R04   R06                                   R18
               y + 3.z + 7.t                    = 7      if you choose  bb = 01               1   3   7                 7                          R05   R07   R09                          R19
                     2.z + 4.t + 6.u           = 1      store these 22 elements                   2   4   6            1          into                    R08   R10   R12                 R20
                              8.t +  u   + 7.v  = 5                                                                 8   1   7       5                                             R11   R13   R15        R21
                                      9.u + 4.v  = 6                                                                      9   4       6                                                       R14   R16        R22

-Then,   1.022  XEQ "3DLS"  >>>>   17.022

  and           R17 = x = -16.47810408            R20 = t = -0.480422463
                  R18 = y =     6.991241632          R21 = u =  0.112313241
                  R19 = z =     1.123905204          R22 = v =  1.247295209

-We have also  R00 = det A = 3882
 

Double Integrals ( equally-spaced arguments )
 

-We want to evaluate  §x1xn  §y1ym  f(x,y).dx.dy    assuming  n & m are odd integers  ( n > 1 , m > 1 )
-The "DIGD" uses Simpson's rule in the directions of  x-axis and y-axis

-  xi  must be  equally spaced:  xi+1 - xi = h     for i = 1,2,...,n-1
-  yj  must be  equally spaced:  yj+1 - yj = k    for j = 1,2,...,m-1
 

Data Registers:       •  R00   = h.k                                          ( Registers R00 thru Rn.m  are to be initialized before executing "DIGD" )

                                  •  R01   = f(x1,y1)  •  Rn+1 = f(x1,y2)  .......  •  Rnm-n+1 = f(x1,ym)
                                  •  R02   = f(x2,y1)  •  Rn+2 = f(x2,y2)  .......  •  Rnm-n+2 = f(x2,ym)
                                     .......................................................................................
                                  •  Rnn   = f(xn,y1)  •  R2n = f(xn,y2)    .......  •  Rnm = f(xn,ym)
Flags:  /
Subroutines:  /
 
 
            STACK             INPUTS          OUTPUTS
                Y                  m                 /
                X                   n  §x1xn §y1ym f(x,y).dx.dy

 
Example:     Calculate   §26 §15  f(x,y).dx.dy   from the following f(x,y) values
 
 
   x\y     1     2     3     4     5
    2     3     4     7     6     3
    4     1     2     4     5     3
    6     4     1     3     4     6

               3  4  7  6  3                   R01  R04  R07  R10  R13
  Store     1  2  4  5  3        in        R02  R05  R08  R11  R14      respectively
               4  1  3  4  6                   R03  R06  R09  R12  R15

-Here, we have  h = 2 and  k = 1   whence  1*2 = 2  STO 00

   n = 3  and  m = 5    whence

     5  ENTER^
     3  XEQ "DIGD"  >>>>   §26 §15  f(x,y).dx.dy  ~  56.8889

-Perfect accuracy is achieved if f is a polynomial and  deg(f) < 4
 

Triple Integrals ( equally-spaced arguments )
 

-Now we evaluate  §x1xn  §y1ym  §z1zp  f(x,y,z).dx.dy.dz    assuming  n , m , p are odd integers  ( n > 1 , m > 1 , p > 1 )
-Simpson's rule in the directions of  x- , y- and z-axis is used.

-  xi  must be  equally spaced:  xi+1 - xi = h         for i = 1,2,...,n-1
-  yj  must be  equally spaced:  yj+1 - yj = h'       for j = 1,2,...,m-1
-  zk  must be  equally spaced:  zk+1 - yk = h''    for j = 1,2,...,p-1
 

Data Registers:       •  R00   = h.h'.h''                                          ( Registers R00 thru Rn.m.p  are to be initialized before executing "DIGT" )

                                  •  R01   = f(x1,y1,z1)  •  Rn+1 = f(x1,y2,z1)  .......  •  Rnm-n+1 = f(x1,ym,z1)
                                  •  R02   = f(x2,y1,z1)  •  Rn+2 = f(x2,y2,z1)  .......  •  Rnm-n+2 = f(x2,ym,z1)
                                     ......................................................................................................
                                  •  Rn     = f(xn,y1,z1)  •  R2n = f(xn,y2,z1)    .......  •  Rnm = f(xn,ym,z1)

                                        •  Rnm+1   = f(x1,y1,z2)  •  Rnm+n+1 = f(x1,y2,z2)  .......  •  R2nm-n+1 = f(x1,ym,z2)
                                        •  Rnm+2   = f(x2,y1,z2)  •  Rnm+n+2 = f(x2,y2,z2)  .......  •  R2nm-n+2 = f(x2,ym,z2)
                                           ......................................................................................................
                                        •  Rnm+n   = f(xn,y1,z2)  •  Rnm+2n = f(xn,y2,z2)    .......  •  R2nm = f(xn,ym,z2)

                                               .............................................................................................................
                                                .............................................................................................................
                                                 .............................................................................................................

                                                   •  Rnm(p-1)+1   = f(x1,y1,zp)  •  Rnm(p-1)+n+1 = f(x1,y2,zp)  .......  •  Rnmp-n+1 = f(x1,ym,zp)
                                                   •  Rnm(p-1)+2   = f(x2,y1,zp)  •  Rnm(p-1)+n+2 = f(x2,y2,zp)  .......  •  Rnmp-n+2 = f(x2,ym,zp)
                                                      ......................................................................................................
                                                   •  Rnm(p-1)+n   = f(xn,y1,zp)  •  Rnm(p-1)+2n = f(xn,y2,zp)    .......  •  Rnmp = f(xn,ym,zp)
Flags:  /
Subroutines:  /
 
 
                STACK                 INPUTS               OUTPUTS
                    Z                      p                      /
                   Y                       m                      /
                    X                      n  §x1xn§y1ym§z1zp f(x,y,z)dx.dy.dz

 
Example:     Evaluate   §13 §15 §17  f(x,y,z).dx.dy.dz   from the following values

-Store:

  f(1,1,1) =  4      f(1,3,1) =  6     f(1,5,1) =  8                    R01   R04   R07
  f(2,1,1) =  7      f(2,3,1) =  9     f(2,5,1) = 11         in        R02   R05   R08     respectively
  f(3,1,1) = 10     f(3,3,1) = 12    f(3,5,1) = 14                    R03   R06   R09

       f(1,1,4) =  64       f(1,3,4) =  96      f(1,5,4) = 128                     R10   R13   R16
       f(2,1,4) = 112      f(2,3,4) = 144     f(2,5,4) = 176           in        R11   R14   R17     respectively
       f(3,1,4) = 160      f(3,3,4) = 192     f(3,5,4) = 224                      R12   R15   R18

            f(1,1,7) =  196      f(1,3,7) = 294     f(1,5,7) = 392                        R19   R22   R25
            f(2,1,7) =  343      f(2,3,7) = 441     f(2,5,7) = 539             in         R20   R23   R26     respectively
            f(3,1,7) =  490      f(3,3,7) = 588     f(3,5,7) = 686                        R21   R24   R27

   h.h'.h'' = 1*2*3 = 6     6  STO 00

  n = m = p = 3   whence   3  ENTER^   ENTER^   XEQ "DIGT"   >>>>   §13 §15 §17  f(x,y,z).dx.dy.dz  =  8208

-These values were actually computed from  f(x,y,z) =  (3x+y).z2   and the result is perfectly accurate since f is a polynomial and deg(f) < 4.
 

Trapezoidal rule ( unequally-spaced arguments too )
 

-We assume a function f is only known by n data points
 
 
       x1         x2    ...........        xn
     f(x1)       f(x2)    ...........      f(xn)

 
-The trapezoidal rule is the easiest formula to evaluate  §x1xn  f(x).dx
-We simply add the areas of n trapezoids:  Sumi=1,2,...,n-1   (xi+1-xi).[ f(xi+1)+f(xi) ]/2
 

Data Registers:       •  R00 = n = number of points ,                        ( Registers R00 thru R2n  are to be initialized before executing "TRAP" )

                                  •  R01 = x1     •  R03 = x2      .......  •  R2n-1 = xn
                                  •  R02 = f(x1)  •  R04 = f(x2)  .......  •  R2n   = f(xn)
Flags:  /
Subroutines:  /
 
 
      STACK        INPUT      OUTPUT
           X              /    §x1xn  f(x).dx

   where n = the number of points.

Example:     Calculate  §18  f(x).dx   from the following values
 
 
      x      1     2.4      4     5.2      7      8
    f(x)      1      4      6      5      4      2

 
-Store these 12 numbers into registers

           R01  R03  R05  R07  R09  R11   ( x )
           R02  R04  R06  R08  R10  R12   ( f(x) )    respectively

-There are 6 points so,     6   STO 00    XEQ "TRAP"    >>>>    §18  f(x).dx  ~  29.2
 

Simpson's rule ( unequally-spaced arguments too )
 

-"SIMP" uses several connected parabolic segments to compute  §x1xn  f(x).dx
-However, if n > 2  is even,  §x1x2  f(x).dx  is calculated by a polynomial of degree 3

-If n = 2 , "TRAP" is called
 

Data Registers:       •  R00 = n = number of points ,                         ( Registers R00 thru R2n  are to be initialized before executing "SIMP" )

                                  •  R01 = x1     •  R03 = x2      .......  •  R2n-1 = xn
                                  •  R02 = f(x1)  •  R04 = f(x2)  .......  •  R2n   = f(xn)
Flags:  /
Subroutines:  /
 
 
      STACK        INPUT      OUTPUT
           X              /    §x1xn  f(x).dx

 where n = the number of points.  ( n > 2 )

Example:     Calculate  §17  f(x).dx  and  §18  f(x).dx   from the following values
 
 
      x      1     2.4      4     5.2      7      8
    f(x)      1      4      6      5      4      2

 
-Store these 12 numbers into registers

    R01  R03  R05  R07  R09  R11   ( x )
    R02  R04  R06  R08  R10  R12   ( f(x) )    respectively

-Five points:   5   STO 00   XEQ "SIMP"  >>>>  §17  f(x).dx  ~  26.4226
-Six points:     6   STO 00         R/S           >>>>  §18  f(x).dx  ~  30.5339
 

Integrals ( connecting cubic segments )
 

-"ITGCUB" uses several segments of cubic curves to evaluate  §x1xn  f(x).dx
  so that the method produces a perfect accuracy if the function is a polynomial of degree < 4

-Registers R00 to R2n are temporarily modified during the calculations, but their contents are finally restored.

-If n = 3 , "SIMP" is used and if n = 2 , "TRAP" is called.
 

Data Registers:       •  R00 = n = number of points ,                          ( Registers R00 thru R2n  are to be initialized before executing "ITGCUB" )

                                  •  R01 = x1     •  R03 = x2      .......  •  R2n-1 = xn
                                  •  R02 = f(x1)  •  R04 = f(x2)  .......  •  R2n   = f(xn)
Flags:  /
Subroutines:  /
 
 
      STACK        INPUT      OUTPUT
           X              /    §x1xn  f(x).dx

  where n = the number of points.  ( n > 3 )

Example:     Calculate  §18  f(x).dx   from the following values
 
 
      x      1     2.4      4     5.2      7      8
    f(x)      1      4      6      5      4      2

 
-Store these 12 numbers into registers

    R01  R03  R05  R07  R09  R11   ( x )
    R02  R04  R06  R08  R10  R12   ( f(x) )    respectively

-There are 6 points,     6   STO 00    XEQ "ITGCUB"    >>>>    §18  f(x).dx  ~  30.2135

Lagrange Interpolation Formula
 

-Given a set of n data points:   M1 ( x1 , y1 )  ,   M2 ( x2 , y2 )  ,  ..........  ,  Mn ( xn , yn )
  store the 2n numbers   x1 , y1 , x2 , y2 , ........ , xn , yn  into contiguous registers,  starting at any register  Rbb:

     •  Rbb = x1       •  Rb+2 = x2      .......     •  Ree-1 = xn
     •  Rb+1 = y1     •  Rb+3 = y2      .......     •  Ree   = yn                  ( ee = bb + 2n -1 )

  LAGR evaluates  L(x)  from x in X-register and the control number of the array bbb.eee in Y-register.
 

REGISTERS:    Rbb thru Ree ,  ee =  2n + bb - 1
FLAGS:            /
 
 
    STACK      INPUTS    OUTPUTS
         Y       bbb.eee      bbb.eee
         X            x        L(x)
         L            /          x

 
Example:

  Evaluate  y(3) and y(5) for the collocation polynomial that takes the values prescribed below:
 
 
    xi     0     1     2     4     7
    yi     3     2     4     6     5

 
-For instance, you store these 10 numbers ( 0 3 1 2 2 4 4 6 7 5 ) into R01 thru R10, then

 1.010  ENTER^
     3     XEQ "LAGR"  produces:  y(3) = 5.847619048

 RDN    5   R/S           ----------   y(5) = 4.523809520

Notes:

-Don't interrupt "LAGR" because synthetic registers P & Q are used
-"LAGR" is called by "LAGRI" hereunder.
 

Integral of the Lagrange Polynomial
 

  "LAGRI"  calculates §x1x2 L(x) dx , given a set of n data points which must be stored as follows:

     •  R00 = n = number of points

     •  R01 = x1     •  R03 = x2      .......  •  R2n-1 = xn
     •  R02 = y1     •  R04 = y2      .......   •  R2n   = yn

REGISTERS:    R00 thru R2n
FLAGS:            /
 
 
      STACK        INPUT      OUTPUT
           X              /    §x1xn  L(x).dx

 
Example:     Calculate  §18  L(x).dx   from the following values
 
 
      x      1     2.4      4     5.2      7      8
      y      1      4      6      5      4      2

 
-Store these 12 numbers into registers   R01  R03  R05  R07  R09  R11   ( x )
                                                            R02  R04  R06  R08  R10  R12    ( y )    respectively
 

-There are 6 points:     6   STO 00      XEQ "LAGRI"     >>>>    §18  L(x).dx  =   29.61789480

Notes:

-Don't interrupt "LAGRI" because synthetic registers P & Q are used
-The data registers are altered by this routine.
 

Natural Cubic Spline Integration
 

  "NCSI"  calculates §x1xn y dx , given a set of n data points which must be stored as follows:
 

REGISTERS:

     •  R00 = n = number of points             R2n+1 to R6n-7  are also used

     •  R01 = x1     •  R03 = x2      .......  •  R2n-1 = xn
     •  R02 = y1     •  R04 = y2      .......   •  R2n   = yn

FLAGS:            /
 
 
      STACK        INPUT      OUTPUT
           X              /    §x1xn  f(x).dx

 
Example:     Calculate  §18  f(x).dx   from the following values
 
 
      x      1     2.4      4     5.2      7      8
    f(x)      1      4      6      5      4      2

 
-Store these 12 numbers into registers   R01  R03  R05  R07  R09  R11   ( x )
                                                            R02  R04  R06  R08  R10  R12   ( f(x) )    respectively

-SIZE 028   ( SIZE = 6n-8 )

-There are 6 points,     6   STO 00    XEQ "NCSI"    >>>>    §18  f(x).dx  =  29.99938860
 

Complete Elliptic Integrals ( 1st & 2nd kinds )
 

-The complete Legendre elliptic integrals are:

 - the complete elliptic integrals of the 1st kind:          K ( m ) =  §0pi/2 ( 1 - m sin2 t )-1/2 dt
 - the complete elliptic integrals of the 2nd kind:         E ( m ) =  §0pi/2  ( 1 - m sin2 t )1/2 dt
 
 
           STACK            INPUTS          OUTPUTS
                Y                 /              E (m)
                X           0 < = m < 1              K (m)

 
Example:

   0.7  XEQ "CEI"  >>>>   E ( 0.7 ) =  2.075363134
                             X<>Y   K ( 0.7 ) = 1.241670567

Note:

-They are computed by the arithmetic-geometric mean ( agm )
 

Incomplete Elliptic Integrals ( 1st & 2nd kinds )
 

-The incomplete Legendre elliptic integrals are:
 

 - the incomplete elliptic integrals of the 1st kind:   F ( v | m ) =  §0v ( 1 - m sin2 t )-1/2 dt       ( 0 < = v < = 90° )
 - the incomplete elliptic integrals of the 2nd kind:  E ( v | m ) =  §0v  ( 1 - m sin2 t )1/2 dt
 
 
           STACK           INPUTS             OUTPUTS
                 T                 /                E ( m ) 
                 Z                 /                 K ( m )
                 Y                m               E ( v | m)
                 X                v ( degrees )               F ( v | m)
                 L                 /               Z ( v | m)

     Z ( v | m)   is in the L-register only if m < 1              ( Z ( v | m )  =  Jacobian Zeta function ).

Examples:

 1-If   v = 84° and m = 0.7

         0.7 ENTER^
         84 XEQ "ELI"  yields       F ( 84° | 0.7 ) = 1.884976271
                                 RDN       E ( 84° | 0.7 ) = 1.184070048
                                 RDN               K ( 0.7 ) = 2.075363134
                                 RDN               E ( 0.7 ) = 1.241670567

               and        LASTX        Z ( 84° | 0.7 ) = 0.056306180

 2-If   v = 84°  and m = 1

            1 ENTER^
           84 XEQ "ELI"  yields     F ( 84° | 1 ) = 2.948700239
                                   RDN      E ( 84° | 1 ) = 0.994521895
                                   RDN              K ( 1 ) = 9.999999999E99    ( in fact infinity )
                                   RDN              E ( 1 ) = 1

Notes:

 K( m ) and E( m ) are automatically calculated in Z- and T registers

 The angle v must be expressed in degrees..
 

Elliptic Integral of the 1st kind
 

-"ELIPF" is an M-Code function written by Angel Martin, which gives the same results as "LEI1"
 

   F ( x | m ) =  §0x ( 1 - m sin2 t ) -1/2 dt  = sin x . RF ( cos2 x ; 1 - m sin2 x ; 1 )
 
 
      STACK        INPUTS      OUTPUT
           Y             m             /
           X             x       F ( x | m )

  >>>   x is expressed in the current angular mode

   in DEG mode

   0.7   ENTER^
   84°  XEQ "ELIPF"  >>>   F( 84° | 0.7 ) = 1.884976271

Note:

  x must be between  -90° ( exclusive ) and +90° ( inclusive )
 

Carlson Elliptic Integral 1st kind , Real Arguments
 

    RF(x,y,z) =  (1/2) §0+infinity  [ ( t + x ).( t + y ).( t + z ) ] -1/2  dt         with  x , y , z  non-negative and at most one is zero
 
 
      STACK        INPUTS      OUTPUTS
           Z             z            /
           Y             y            /
           X             x      RF(x,y,z)
           L             /            x

 
    2   ENTER^
    3   ENTER^
    4   XEQ "RF"  >>>   RF(4,3,2) = 0.584082842
 

Carlson Elliptic Integral 1st kind , Conjugate Complex Arguments
 

  "RFZ"  calculates  RF(x,y+i.z,y-i.z)
 
 
      STACK        INPUTS      OUTPUTS
           Z             z            /
           Y             y            /
           X             x   RF(x,y+i.z,y-i.z)
           L             /            x

 
  2   ENTER^
  3   ENTER^
  4   XEQ "RFZ"  >>>   RF(4,3+2.i,3-2.i) = 0.529597761
 

Carlson Elliptic Integral 3rd kind , Real Arguments
 

    RJ(x,y,z,p) =  (3/2)  §0+infinity  ( t + p ) -1 [ ( t + x ).( t + y ).( t + z ) ] -1/2   dt        with  x , y , z  non-negative and at most one is zero
 
 
      STACK        INPUT      OUTPUT
           T          p > 0            /
           Z             z            /
           Y             y            /
           X             x     RJ(x,y,z,p)

 
  1   ENTER^
  2   ENTER^
  3   ENTER^
  4   XEQ "RJ"  >>>   RJ(4,3,2,1) = 0.360378094

-Alpha register is cleared.
 

Carlson Elliptic Integral 3rd kind , Conjugate Complex Arguments
 

    "RJZ" computes  RJ(x,y+i.z,y-i.z,p)  with p > 0
 
 
      STACK        INPUTS      OUTPUT
           T          p > 0            /
           Z             z            /
           Y             y            /
           X             x  RJ(x,y+i.z,y-i.z,p)

 
  1   ENTER^
  2   ENTER^
  3   ENTER^
  4   XEQ "RJZ"  >>>   RJ(4,3+2.i,3-2.i,1) = 0.287311651

-Alpha register is cleared.
 

Carlson Symmetric Elliptic Integral 2nd kind , Real Arguments
 

     RG(x,y,z) =  (1/4)  §0+infinity  [  ( t + x ).( t + y ).( t + z ) ] -1/2 [ x/(t+x) + y/(t+y) + z/(t+z) ]  t.dt
 

REGISTERS:    R00 thru R03
FLAGS:            /
 
 
      STACK        INPUTS      OUTPUT
           Z             z            /
           Y             y            /
           X             x      RG(x,y,z)

 
  2   ENTER^
  3   ENTER^
  4   XEQ "RG"  >>>   RG(4,3,2) = 1.725503029
 

Legendre Elliptic Integral of the 1st kind
 

    F ( x | m ) =  §0x ( 1 - m sin2 t ) -1/2 dt  = sin x . RF ( cos2 x ; 1 - m sin2 x ; 1 )
 

REGISTERS:    R00 thru R03
FLAGS:            F01-F02
 
 
      STACK        INPUTS      OUTPUT
           Y             m             /
           X             x       F ( x | m )

  >>>   x is expressed in the current angular mode

   in DEG mode

   0.7   ENTER^
   84°  XEQ "LEI1"  >>>   F( 84° | 0.7 ) = 1.884976271

Note:

  x must be between  -90° ( exclusive ) and +90° ( inclusive )
 

Legendre Elliptic Integral of the 2nd kind
 

     E ( x | m ) =  §0x ( 1 - m sin2 t ) 1/2 dt  = sin x . RF ( cos2 x ; 1 - m.sin2 x ; 1 ) - (m/3) sin3 x . RJ ( cos2 x ; 1 - m.sin2 x ; 1 ; 1 )
 

REGISTERS:    R00 thru R04
FLAGS:            F01-F02
 
 
      STACK        INPUTS      OUTPUT
           Y             m             /
           X             x       E ( x | m )

  >>>   x is expressed in the current angular mode

   in DEG mode

   0.7   ENTER^
   84°  XEQ "LEI2"  >>>   1.184070048

Note:

  x must be between  -90° ( exclusive ) and +90° ( inclusive )
 

Legendre Elliptic Integral of the 3rd kind
 

     ¶ ( n ; x | m ) =  §0x  ( 1 + n sin2 t ) -1 ( 1 - m sin2 t ) -1/2 dt
                         = sin x . RF ( cos2 x ; 1 - m.sin2 x ; 1 ) - (n/3) sin3 x . RJ ( cos2 x ; 1 - m.sin2 x ; 1 ; 1 + n.sin2 x )
 

REGISTERS:    R00 thru R04
FLAGS:            F01-F02
 
 
      STACK        INPUTS      OUTPUT
           Z             n             /
           Y             m             /
           X             x    ¶ ( n ; x | m )

  >>>   x is expressed in the current angular mode

   in DEG mode

   0.9   ENTER^
   0.7   ENTER^
   84°  XEQ "LEI3"  >>>   1.336853616

Note:

  x must be between  -90° ( exclusive ) and +90° ( inclusive )
-The sign convention for n is opposite that of Abramowitz & Stegun
 

Arc length of a curve
 

-The arc length of the curve  y = f(x)   ( a < x < b )  is given by   L = §ab  ( 1 + (y')2 )1/2 dx
-The Romberg extrapolation to the limit is used to compute this integral without calculating the derivative.
 

Data Registers:

     R00:  f name
     R01:  a
     R02:  b

Flag:   F02
Subroutines:   "ROM" and a program which takes x in X-register and calculates f(x) in X-register.
 
 
      STACK        INPUTS      OUTPUTS
        alpha         f name            /
           Y             a            /
           X             b      Arc Length

 
Example:

-Calculate the arc length of the curve   y = ln x      1 < x < 3
 
 
    LBL "T"
    LN
    RTN

 
  alpha   "T"  alpha     FIX 9

  1    ENTER^
  3    XEQ   "CURVE"     the following approximations are displayed:

    2.300459681
    2.301931038
    2.301986835
    2.301987537
    2.301987533

-The exact result is  2.301987535  ( rounded to 9 decimals )
 

Romberg Integration - midpoint formula - Step Doubling
 

   I  =   §a f(x) dx
 

-The following program uses the very simple "midpoint" formula:   §ab f(x) dx = (b-a).f((a+b)/2)
 

Data Registers:

     R00:  f name
     R01:  a
     R02:  b

Flags: /
Subroutines:  "ROM" and a program which takes x in X-register and calculates f(x) in X-register.
 
 
 
      STACK        INPUTS      OUTPUTS
        alpha         f name             /
           Y             a             /
           X             b     §a f(x) dx

 
Example:

-Evaluate  I =  §12  xx  dx   to  7 decimals
 
 
  LBL "T"
  ENTER^
  Y^X
  RTN

 
alpha   "T"   alpha     FIX 7

   1   ENTER^
   2   XEQ  "IG"    the following values are displayed:

  2.0438805
  2.0503885
  2.0504461
  2.0504462

-In fact, the result is correct to almost 9 decimals:  I = 2.050446234  ( the last digit should be a 5 )
 

Romberg Integration - midpoint formula - Step Tripling
 

   I  =   §a f(x) dx
 

-Like "IG" , "IG3"  uses the midpoint formula,
 but in order to use the previous evaluations of the function,
 the number of steps is trippled with each iteration.
 

Data Registers:            R00 thru R09 are unused

      R10:  f name       R13:  3n                  R16:  counter and 32 , 92 , 272 ,  ....
      R11:  a                R14:  sums              R17:  counter
      R12:  (b-a)/3n     R15:  x                    R18, R19 , .... = successive approximations

Flags: /
Subroutine:   a program which calculates f(x) assuming x is in X-register upon entry.
 
 
 
      STACK        INPUTS      OUTPUTS
        alpha         f name             /
           Y             a             /
           X             b      §a f(x) dx

 
Example:    Evaluate   §13  x ( 1+x3 )1/2  dx
 
 
  LBL "T"
  ENTER^
  ENTER^
  X^2
  *
  ENTER^
  SIGN
  +
  SQRT
  *
  RTN

 
-The accuracy is controlled by the display setting

  alpha   "T"   alpha    FIX 7

         1   ENTER^
         3   XEQ "IG3"   >>>>   13.7718432
                                              13.7693525
                                              13.7693320
                                              13.7693320

-The last value is  13.76933201  ( the last decimal should be a "2" )
-Actually, the next to last result was exact:  13.76933202  ( in R15 and in L-register )

-Unlike "IG" , "IG3" doesn't lose the previous function evaluations, but step trippling often leads to long execution times.
 

Midpoint formula - Double Integrals
 

   I  =   §ab §u(x)v(x)  f(x;y) dx dy

-The mid-point formula is used in both x and y-axis.
 

Data Registers:                                             ( Registers R04-R05 are to be initialized before executing "IGD" )

     R00:  f name        R03:  n
     R01:  a             •  R04:    u name
     R02:  b             •  R05:    v name

Flags: /
Subroutines:   "ROM"  and

   A program which takes x in X-register and y in Y-register and calculates f(x,y) in X-register.
   A program which takes x in X-register and calculates u(x) in X-register.
   A program which takes x in X-register and calculates v(x) in X-register.
 
 
      STACK        INPUTS      OUTPUTS
        alpha         f name             /
           Y             a             /
           X             b        Integral

 
Example:

-Calculate    §23 §xx^2  ( 1 + xy )1/2 dx dy   to  7 decimals
 
 
  LBL "T"
  *
  1
  +
  SQRT
  RTN
  LBL "U"
  RTN
  LBL "V"
  X^2
  RTN

 
  "U"     ASTO 04
  "V"     ASTO 05     FIX 7

   alpha  "T"  alpha

  2   ENTER^
  3   XEQ  "IGD"    the HP-41 displays:

  13.7800635
  13.7747491
  13.7746576
  13.7746565
 

Midpoint formula - Triple Integrals
 

  I  =  §ab §u(x)v(x)  §w(x,y)t(x,y)  f(x,y,z) dx dy dz
 

-The mid-point formula is again used ( in x, y and z-axis ).
 

Data Registers:                                       ( Registers R04-R05-R06-R07 are to be initialized before executing "IGT" )

         R00:  f name     •  R04:    u name
         R01:  a             •  R05:    v name
         R02:  b             •  R06:   w name
         R03:  n             •  R07    t name

Flags: /
Subroutines:   "ROM"

  A program which takes x in X-register, y in Y-register and z in Z-register and calculates f(x,y,z) in X-register.
  A program which takes x in X-register and calculates u(x) in X-register.
  A program which takes x in X-register and calculates v(x) in X-register.
  A program which takes x in X-register and y in Y-register and calculates w(x,y) in X-register.
  A program which takes x in X-register and y in Y-register and calculates t(x,y) in X-register.
 
 
      STACK        INPUTS      OUTPUTS
        alpha         f name             /
           Y             a             /
           X             b       Integral

 
Example:   Evaluate    §01 §0x §-x-y-y   1 / ( 1 + x + y + z )   dx dy dz   to 5 places
 
 
  LBL "P"
  +
  +
  1
  +
  1/X
  RTN
  LBL "U"
  CLX
  RTN
  LBL "V"
  RTN
  LBL "W"
  +
  CHS
  RTN
  LBL "T"
  X<>Y
  CHS
  RTN

 
   "U"     ASTO 04     "V" ASTO 05
   "W"    ASTO 06     "T"  ASTO 07

  alpha  "P"  alpha   FIX 5

    0  ENTER^
    1  XEQ  "IGT"    the calculator displays:

  0.24838
  0.24998
  0.25000                 ( exact value = 1/4 )

-The last approximation is actually    0.249999724
 

A subroutine: "ROM"
 

-This subroutine is called by many programs in the Romberg section.
-The calling program must have the following structure:

  LBL "???"     LBL 01       calculates the               XEQ "ROM"     RCL IND 22
  1                   ..........            integral                     X#0?                 END
  STO 03        ...........        corresponding to n       GTO 01

-The accuracy is controlled by the display setting.

-Registers R21-R22- .... are used and  R03 = n
 

Area of a Surface of Revolution
 

-The rotation of the curve  y = f(x)   ( a < x < b )  around x-axis generates a surface of revolution given by   S = 2.pi  §ab  y ( 1 + (y')2 )1/2 dx
  "SREVL" doesn't use this formula and avoids the calculation of dy/dx : the area of a truncated cone is used instead.
 

Data Registers:

     R00:  f name
     R01:  a
     R02:  b

Flag:   F02
Subroutines:  "ROM" and a program which takes x in X-register and calculates f(x) in X-register.
 
 
      STACK        INPUTS      OUTPUTS
        alpha         f name             /
           Y             a             /
           X             b         Area

 
Example:

-Evaluate the area of the surface of revolution generated by the rotation of the curve   y = sin x  ( 0 < x < pi )  around the x-axis.
 
 
  LBL "T"
  SIN
  RTN

 
  XEQ RAD    FIX 9

  alpha  "T"  alpha

   0    ENTER^
  PI   XEQ   "SREVL"     the following approximations are displayed:

  15.59985804
  14.26493243
  14.42990383
  14.42356623
  14.42359884
  14.42359950
 

PPC ROM Integrator
 

-"PPCIG" is the program given in the PPC ROM. It employs a Romberg method to evaluate      I  =   §a f(x) dx

1°)  Place the name of the subroutine in the alpha register ( the function name )  and  ASTO 10
2°)  Place a ENTER^ b  ( i-e  a in Y-register and b in X-register )
3°)  XEQ "PPCIG"

-The accuracy is controlled by the display setting.
-If flag F10 is set, the HP-41 will display the successive approximations.
 
 
      STACK        INPUTS      OUTPUTS
           Y             a             /
           X             b      §ab  f(x)  dx 

 
Another Integrator
 

-"ITG" also computes  §a f(x) dx
 
 
      STACK        INPUTS      OUTPUTS
        Alpha         F name             /
           Z             a             /
           Y             b             /
           X             n      §ab  f(x)  dx 

 
Example:     Evaluate  §13  e-x^2 dx

-Load this short routine:
 
 
 01  LBL "T" 
 02  X^2
 03  CHS
 04  E^X
 05  RTN

 
  Alpha  "T"  Alpha

  1    ENTER^
  3    ENTER^                                                    and with  n = 8
  8    XEQ "ITG"  >>>>  0.139381955
 

Integral from a to b
 

-"IGAB"  evaluates  I  =   §a f(x) dx  by Gauss-Legendre 16-point formula over  [a,b] without subdivising the interval
-This program prompts for each input ( fname , a , b )
 

Integral from a to infinity
 

-"IGAI"  evaluates  I  =   §a+infinity  f(x) dx  by Gauss-Legendre 16-point formula after a change of variable to change the interval into a bounded one.
-This program prompts for each input ( fname , a )
 

PPC Solver
 

-"PPCSV" solves non-linear equations of the type  f(x) = 0   ( it comes from the PPC ROM )
-It employs the secant method and the accuracy is controlled by the display format.
 
 
      STACK        INPUTS      OUTPUTS
        Alpha        F name             /
           Y             a             /
           X             b         Root 

 
Example:    Find a solution of  x3 - 3 x + 1 = 0   near   1 & 2
 
 
  LBL "W"
  ENTER^
  X^2
  3
  -
  *
  1
  +
  RTN

 
  Alpha  "W"  Alpha       FIX 9

  1   ENTER^
  2   XEQ "PPCSV"   >>>>   Root = 1.532088886
 

Differential Equations
 

 "DIFEQ"  solves 1st and 2nd order differential equations
 

Example:     y" + 2 x y' + y = 0   i-e   y" = - 2 x y' - y  with initial values  y(0) = 1 ,  y'(0) = 0

->  Compute  y(0.1)
 
 
  LBL "T"
  RCL Z
  *
  ST+ X
  +
  CHS
  RTN

 
   XEQ "DIFEQ"   >>>>    " NAME?"

    T     R/S    >>>>    "ORDER=?"
    2     R/S    >>>>    "STEP SIZE=?"
  0.1    R/S    >>>>    "X0=?"
   0      R/S    >>>>    "Y0=?"
   1      R/S    >>>>    "Y0.=?"
   0      R/S    >>>>     0.1 = R02                R/S         0.995021 = y(0.1) = R03   and we have  y' (0.1) = -0.099170 = R04

-Press  R/S  to  get  y(0.2)  &  y' (0.2)   .... and so on ...
 

Data Points for Discrete Integration
 

-XEQ "RGDATA" to store the data points for discrete cases.
 
 

References:

[1]   "Numerical Analysis" - F. Scheid - Mc Graw Hill  ISBN 07-055197-9
[2]   "PPC-ROM user's manual"  ( page 220 ..... )
[3]   "Extended Functions made easy" - Keith Jarett ( synthetix )
[4]   "Numerical Recipes in C++" - Press , Vetterling , Teukolsky , Flannery  - Cambridge University Press - ISBN  0-521-75033-4