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 |
-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 |
-"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
-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 = §ab
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 | §ab 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 = §ab
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 | §ab 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
-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
-"PPCIG" is the program given in the PPC ROM. It employs a Romberg method to evaluate I = §ab 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 |
-"ITG" also computes §ab 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
-"IGAB" evaluates I = §ab
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 )
-"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 )
-"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
"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