hp41programs

INTEGRALEQ

Integral Equations for the HP-41


Overview
 

 1°)  Fredholm Equation of the 2nd kind over [a,b]

  a) With Gauss-Legendre 10-point formula
  b) With Gauss-Legendre 15-point formula

 2°)  Fredholm Equation of the 2nd kind over [0,+infinity[

  a) With Gauss-Laguerre 10-point formula
  b) With Gauss-Laguerre 15-point formula

 3°)  Fredholm Equation of the 2nd kind over ]-infinity, +infinity[

  a) With Gauss-Hermite 10-point formula
  b) With Gauss-Hermite 15-point formula

 4°)  Volterra Equation of the 1st kind
 5°)  Volterra Equation of the 2nd kind
 

-Fredholm equations of the 2nd kind are solved by the Nystrom method:
-We solve a linear system that is usually well-conditionned.
-With Fredholm equations of the 1st kind, the same method often leads to badly ill-conditionned systems, so they are omitted here.
 

1°) Fredholm Integral Equation 2nd kind over [a,b]
 

  a) With Gauss-Legendre 10-point formula
 

 "2FHAB" solves the integral equation:

   f(x) = g(x) + §ab K(x,y) f(y) dy    where a , b , g & K are known and f is the unknown function.

-Gauss-Legendre 10-point formula leads to a linear system of 10 equations in 10 unknowns

-After solving this system, the integral equation is used again to evaluate f(x) for a given x
-With an interpolation formula, we would lose the accuracy given by Gaussian integration.
 

Data Registers:              R00 = det              ( Registers R01 & R02 are to be initialized before executing "2FHAB" )

                                      •  R01 = a                  R03 = x                   R11 to R30: abscissas & weights
                                      •  R02 = b                  R04 to R06: temp    R31 to R140: coeff of the system

Flags: /
Subroutines:  "LS"  ( cf "Linear & Non-Linera Systems for the HP-41" )
                         "K"   which must take x & y in registers X & Y and return K(x,y) in X-register
                        "GG" ------------------  in X-register and return g(x) in X-register
 
 

 01 LBL "2FHAB"
 02 STO 03
 03 130
 04 STO 04
 05 10
 06 STO 05
 07 STO 06
 08 .148874339
 09 RCL 02
 10 RCL 01              
 11 -
 12 2
 13 /
 14 *
 15 STO 11
 16 CHS
 17 STO 16
 18 .4333953941
 19 LASTX
 20 *
 21 STO 12
 22 CHS
 23 STO 17
 24 .6794095683
 25 LASTX
 26 *
 27 STO 13
 28 CHS
 29 STO 18
 30 .8650633667
 31 LASTX
 32 *
 33 STO 14
 34 CHS
 35 STO 19
 36 .9739065285
 37 LASTX
 38 *
 39 STO 15              
 40 CHS
 41 STO 20
 42 .2955242247
 43 LASTX
 44 *
 45 STO 21
 46 STO 26
 47 .2692667193
 48 LASTX
 49 *
 50 STO 22
 51 STO 27
 52 .2190863625
 53 LASTX
 54 *
 55 STO 23
 56 STO 28
 57 .1494513492
 58 LASTX
 59 *
 60 STO 24
 61 STO 29
62 6.667134431 E-2
 63 LASTX
 64 *
 65 STO 25              
 66 STO 30
 67 RCL 01
 68 RCL 02
 69 +
 70 2
 71 /
 72 ST+ 11
 73 ST+ 12
 74 ST+ 13
 75 ST+ 14
 76 ST+ 15
 77 ST+ 16
 78 ST+ 17
 79 ST+ 18
 80 ST+ 19
 81 ST+ 20
 82 LBL 01
 83 RCL 06
 84 RCL 05
 85 10
 86 ST+ Z
 87 +
 88 RCL IND Y
 89 RCL IND Y
 90 XEQ "K"
 91 RCL 06
 92 20
 93 +
 94 X<>Y
 95 RCL IND Y
 96 *
 97 RCL 05              
 98 RCL 06
 99 X=Y?
100 GTO 00
101 X<> Z
102 GTO 01
103 LBL 00
104 RDN
105 SIGN
106 -
107 LBL 01
108 STO IND 04
109 DSE 04
110 DSE 05
111 GTO 01
112 10
113 STO 05
114 DSE 06
115 GTO 01
116 140
117 STO 04
118 LBL 02
119 RCL 05
120 10
121 +
122 RCL IND X
123 XEQ "GG"
124 CHS
125 STO IND 04
126 DSE 04
127 DSE 05
128 GTO 02
129 31.1401
130 ENTER
131 CLX
132 XEQ "LS"
133 1/X
134 RCL 03             
135 LBL 10
136 STO 03
137 140
138 STO 04
139 20.01
140 STO 05
141 CLX
142 STO 06
143 LBL 03
144 RCL IND 05
145 RCL 03
146 XEQ "K"
147 RCL 05
148 10
149 +
150 X<>Y
151 RCL IND Y
152 *
153 RCL IND 04
154 *
155 ST+ 06
156 DSE 04
157 DSE 05
158 GTO 03
159 RCL 03
160 XEQ "GG"
161 ST+ 06
162 RCL 03             
163 SIGN
164 RCL 06
165 RTN
166 GTO 10
167 END

 
    ( 366 bytes / SIZE 141 )
 
 

      STACK        INPUTS      OUTPUTS
           X             x           f(x)
           L             /            x

 
Example:     Find f such that  f(x) = exp(-x^2) + (2/e) Sinh x + §-1+1 ( -x - 2.y ) exp(-x.y) f(y) dy

-Load these 2 subroutines in main memory
 
 

 01 LBL "GG"
 02 ENTER
 03 E^X-1
 04 LASTX
 05 CHS
 06 E^X-1
 07 -
 08 1
 09 E^X
 10 /
 11 X<>Y
 12 X^2
 13 CHS
 14 E^X
 15 +
 16 RTN
 17 LBL "K"
 18 RCL Y
 19 RCL Y
 20 *
 21 CHS
 22 E^X
 23 X<>Y
 24 R^
 25 ST+ X
 26 +
 27 *
 28 CHS
 29 RTN
 30 END

 
  1  STO 02  CHS  STO 01  and if you want to calculate f(0.7)

  0.7  XEQ "2FHAB"  >>>>  f(0.7) = 0.612626396                               ---Execution time = 7m21s---

-After executing "2FHAB", if you want to calculate another f(x), place the new x-value in X-register and press R/S or XEQ 10

   1  R/S  >>>>   f(1) = 0.367879441                                   ---Execution time = 14s---
   0  R/S  >>>>   f(0) = 1.000000005

   3  R/S  >>>>   f(3) = 0.000123336  ( exact result = 0.000123410... )

Notes:

-You could also press RTN R/S or XEQ "2FHAB" again, but it would re-solve the linear system which obviously takes a long time.

-Since the exact solution is f(x) = exp(-x^2) , the precision is very good... provided x remains in the interval [a,b]
 

  b) With Gauss-Legendre 15-point formula
 

 "2FHAB+" solves the same integral equation:

   f(x) = g(x) + §ab K(x,y) f(y) dy    where a , b , g & K are known and f is the unknown function.

-But Gauss-Legendre 15-point formula is used which leads to a linear system of 15 equations in 15 unknowns

-After solving this system, the integral equation is used again to evaluate f(x) for a given x
 

Data Registers:              R00 = det             ( Registers R01 & R02 are to be initialized before executing "2FHAB+" )

                                      •  R01 = a                  R03 = x                   R07 to R36: abscissas & weights
                                      •  R02 = b                  R04 to R06: temp    R37 to R276: coeff of the system

Flags: /
Subroutines:  "LS"  ( cf "Linear & Non-Linera Systems for the HP-41" )
                         "K"   which must take x & y in registers X & Y and return K(x,y) in X-register
                        "GG" ------------------  in X-register and return g(x) in X-register
 
 

 01 LBL "2FHAB+"
 02 STO 03
 03 261
 04 STO 04
 05 15
 06 STO 05
 07 STO 06
 08 CLX
 09 STO 07
 10 .201194094
 11 RCL 02
 12 RCL 01              
 13 -
 14 2
 15 /
 16 *
 17 STO 08
 18 CHS
 19 STO 09
 20 .3941513471
 21 LASTX
 22 *
 23 STO 10
 24 CHS
 25 STO 11
 26 .5709721726
 27 LASTX
 28 *
 29 STO 12
 30 CHS
 31 STO 13
 32 .7244177314
 33 LASTX
 34 *
 35 STO 14
 36 CHS
 37 STO 15
 38 .8482065834
 39 LASTX
 40 *
 41 STO 16
 42 CHS
 43 STO 17
 44 .9372733924
 45 LASTX
 46 *
 47 STO 18
 48 CHS
 49 STO 19              
 50 .987992518
 51 LASTX
 52 *
 53 STO 20
 54 CHS
 55 STO 21
 56 .2025782419
 57 LASTX
 58 *
 59 STO 22
 60 .1984314853
 61 LASTX
 62 *
 63 STO 23
 64 STO 24
 65 .186161
 66 LASTX
 67 *
 68 STO 25
 69 STO 26
 70 .1662692058
 71 LASTX
 72 *
 73 STO 27
 74 STO 28
 75 .1395706779
 76 LASTX
 77 *
 78 STO 29
 79 STO 30
 80 .1071592205
 81 LASTX
 82 *
 83 STO 31
 84 STO 32              
 85 7.036604749 E-2
 86 LASTX
 87 *
 88 STO 33
 89 STO 34
 90 3.0753242 E-2
 91 LASTX
 92 *
 93 STO 35
 94 STO 36
 95 RCL 01
 96 RCL 02
 97 +
 98 2
 99 /
100 ST+ 07
101 ST+ 08
102 ST+ 09
103 ST+ 10
104 ST+ 11
105 ST+ 12
106 ST+ 13
107 ST+ 14
108 ST+ 15
109 ST+ 16
110 ST+ 17
111 ST+ 18
112 ST+ 19
113 ST+ 20
114 ST+ 21
115 LBL 01
116 RCL 06              
117 RCL 05
118 6
119 ST+ Z
120 +
121 RCL IND Y
122 RCL IND Y
123 XEQ "K"
124 RCL 06
125 21
126 +
127 X<>Y
128 RCL IND Y
129 *
130 RCL 05
131 RCL 06
132 X=Y?
133 GTO 00
134 X<> Z
135 GTO 01
136 LBL 00
137 RDN
138 SIGN
139 -
140 LBL 01
141 STO IND 04
142 DSE 04
143 DSE 05
144 GTO 01
145 15
146 STO 05
147 DSE 06
148 GTO 01
149 276
150 STO 04             
151 LBL 02
152 RCL 05
153 6
154 +
155 RCL IND X
156 XEQ "GG"
157 CHS
158 STO IND 04
159 DSE 04
160 DSE 05
161 GTO 02
162 37.27615
163 ENTER
164 CLX
165 XEQ "LS"
166 1/X
167 RCL 03
168 LBL 10
169 STO 03
170 276
171 STO 04
172 21.006
173 STO 05             
174 CLX
175 STO 06
176 LBL 03
177 RCL IND 05
178 RCL 03
179 XEQ "K"
180 RCL 05
181 15
182 +
183 X<>Y
184 RCL IND Y
185 *
186 RCL IND 04
187 *
188 ST+ 06
189 DSE 04
190 DSE 05
191 GTO 03
192 RCL 03
193 XEQ "GG"
194 ST+ 06
195 RCL 03
196 SIGN
197 RCL 06
198 RTN
199 GTO 10
200 END

 
    ( 457 bytes / SIZE 277 )
 
 

      STACK        INPUTS      OUTPUTS
           X             x           f(x)
           L             /            x

 
Example:  The same one:   find f such that  f(x) = exp(-x^2) + (2/e) Sinh x + §-1+1 ( -x - 2.y ) exp(-x.y) f(y) dy

-Load these 2 subroutines in main memory
 
 

 01 LBL "GG"
 02 ENTER
 03 E^X-1
 04 LASTX
 05 CHS
 06 E^X-1
 07 -
 08 1
 09 E^X
 10 /
 11 X<>Y
 12 X^2
 13 CHS
 14 E^X
 15 +
 16 RTN
 17 LBL "K"
 18 RCL Y
 19 RCL Y
 20 *
 21 CHS
 22 E^X
 23 X<>Y
 24 R^
 25 ST+ X
 26 +
 27 *
 28 CHS
 29 RTN
 30 END

 
  1  STO 02  CHS  STO 01  and if you want to calculate f(0.7)

  0.7  XEQ "2FHAB+"  >>>>  f(0.7) = 0.612626397                               ---Execution time = 18m54s---

-After executing "2FHAB+", if you want to calculate another f(x), place the new x-value in X-register and press R/S or XEQ 10

   1  R/S  >>>>   f(1) = 0.367879440                                   ---Execution time = 20s---
   0  R/S  >>>>   f(0) = 1.000000008

   3  R/S  >>>>   f(3) = 0.000123287  ( exact result = 0.000123410... )

Notes:

-You could also press RTN R/S or XEQ "2FHAB+" again, but it would re-solve the linear system which obviously takes a long time.

-Since the exact solution is f(x) = exp(-x^2) , the precision is very good... provided x remains in the interval [a,b],
 though not better than with "2FHAB" in this example.
 

2°)  Fredholm Integral Equation 2nd kind over [0,+infinity[
 

  a) With Gauss-Laguerre 10-point formula
 

 "2FH0I" solves the integral equation:

   f(x) = g(x) + §0+infinity K(x,y) f(y) dy    where a , b , g & K are known and f is the unknown function.

-We could make a change of argument ( like y = Tan z ) and use "2FHAB"
-But "2FH0I" employs Gauss-Laguerre 10-point formula which leads to a linear system of 10 equations in 10 unknowns

-After solving this system, the integral equation is used again to evaluate f(x) for a given x
 

Data Registers:              R00 = det                R01-R02: unused

                                         R03 = x                   R11 to R30: abscissas & weights
                                         R04 to R06: temp    R31 to R140: coeff of the system

Flags: /
Subroutines:  "LS"  ( cf "Linear & Non-Linera Systems for the HP-41" )
                         "K"   which must take x & y in registers X & Y and return K(x,y) in X-register
                        "GG" ------------------  in X-register and return g(x) in X-register
 
 

 01 LBL "2FH0I"
 02 STO 03
 03 130
 04 STO 04
 05 10
 06 STO 05
 07 STO 06
 08 .1377934705
 09 STO 11
 10 .7294545495
 11 STO 12
 12 1.808342902
 13 STO 13
 14 3.401433698
 15 STO 14
 16 5.55249614
 17 STO 15
 18 8.330152747
 19 STO 16
 20 11.84378584
 21 STO 17
 22 16.27925783
 23 STO 18
 24 21.99658581
 25 STO 19
 26 29.92069701
 27 STO 20
 28 .3540097386
 29 STO 21
 30 .831902301
 31 STO 22
 32 1.330288562
 33 STO 23
 34 1.863063903
 35 STO 24
 36 2.450255558
 37 STO 25        
 38 3.122764155
 39 STO 26
 40 3.934152696
 41 STO 27
 42 4.992414872
 43 STO 28
 44 6.572202485
 45 STO 29
 46 9.78469584
 47 STO 30
 48 LBL 01
 49 RCL 06
 50 RCL 05
 51 10
 52 ST+ Z
 53 +
 54 RCL IND Y
 55 RCL IND Y
 56 XEQ "K"
 57 RCL 06        
 58 20
 59 +
 60 X<>Y
 61 RCL IND Y
 62 *
 63 RCL 05
 64 RCL 06
 65 X=Y?
 66 GTO 00
 67 X<> Z
 68 GTO 01
 69 LBL 00
 70 RDN
 71 SIGN
 72 -
 73 LBL 01
 74 STO IND 04
 75 DSE 04
 76 DSE 05
 77 GTO 01
 78 10
 79 STO 05        
 80 DSE 06
 81 GTO 01
 82 140
 83 STO 04
 84 LBL 02
 85 RCL 05
 86 10
 87 +
 88 RCL IND X
 89 XEQ "GG"
 90 CHS
 91 STO IND 04
 92 DSE 04
 93 DSE 05
 94 GTO 02
 95 31.1401
 96 ENTER
 97 CLX
 98 XEQ "LS"
 99 1/X
100 RCL 03
101 LBL 10         
102 STO 03
103 140
104 STO 04
105 20.01
106 STO 05
107 CLX
108 STO 06
109 LBL 03
110 RCL IND 05
111 RCL 03
112 XEQ "K"
113 RCL 05
114 10
115 +
116 X<>Y
117 RCL IND Y
118 *
119 RCL IND 04
120 *
121 ST+ 06
122 DSE 04
123 DSE 05
124 GTO 03
125 RCL 03        
126 XEQ "GG"
127 ST+ 06
128 RCL 03
129 SIGN
130 RCL 06
131 RTN
132 GTO 10
133 END

 
    ( 417 bytes / SIZE 141 )
 
 

      STACK        INPUTS      OUTPUTS
           X             x           f(x)
           L             /            x

 
Example:    Find f such that  f(x) = exp(-x) - x + §0+infinity x y f(y) dy
 
 

 01 LBL "GG"
 02 CHS
 03 E^X
 04 LASTX
 05 +
 06 RTN
 07 LBL "K"
 08 *
 09 RTN
 10 END

 
-If you want to evaluate f(0.7)

  0.7  XEQ "2FH0I"  >>>>  0.496585304                                          ---Execution time = 6m40s---

-Likewise,

      0  R/S ( or XEQ 10 )  >>>>  1.000000000                        ---Execution time = 9s---
      1  R/S ( or XEQ 10 )  >>>>  0.367879442

-The precision is again very good.
 

  b) With Gauss-Laguerre 15-point formula
 

 "2FH0I+" solves the same integral equation,   f(x) = g(x) + §0+infinity K(x,y) f(y) dy    but with Gauss-Laguerre 15-point formula:
 

Data Registers:              R00 = det           R01-R02: unused

                                         R03 = x                   R07 to R36: abscissas & weights
                                         R04 to R06: temp    R37 to R276: coeff of the system

Flags: /
Subroutines:  "LS"  ( cf "Linear & Non-Linera Systems for the HP-41" )
                         "K"   which must take x & y in registers X & Y and return K(x,y) in X-register
                        "GG" ------------------  in X-register and return g(x) in X-register
 
 

 01 LBL "2FH0I+"
 02 STO 03
 03 261
 04 STO 04
 05 15
 06 STO 05
 07 STO 06
 08 9.330781202 E-2
 09 STO 07              
 10 .4926917403
 11 STO 08
 12 1.215595412
 13 STO 09
 14 2.269949526
 15 STO 10
 16 3.667622722
 17 STO 11
 18 5.425336627
 19 STO 12
 20 7.565916227
 21 STO 13
 22 10.12022857
 23 STO 14
 24 13.13028248
 25 STO 15
 26 16.65440771
 27 STO 16
 28 20.7764789
 29 STO 17
 30 25.62389423
 31 STO 18
 32 31.40751917
 33 STO 19
 34 38.53068331
 35 STO 20              
 36 48.02608557
 37 STO 21
 38 .2395781703
 39 STO 22
 40 .5601008428
 41 STO 23
 42 .8870082629
 43 STO 24
 44 1.223664402
 45 STO 25
 46 1.574448722
 47 STO 26
 48 1.944751977
 49 STO 27
 50 2.341502057
 51 STO 28
 52 2.774041927
 53 STO 29
 54 3.255643346
 55 STO 30
 56 3.806311714
 57 STO 31
 58 4.458477754
 59 STO 32              
 60 5.270017784
 61 STO 33
 62 6.35956347
 63 STO 34
 64 8.031787632
 65 STO 35
 66 11.5277721
 67 STO 36
 68 LBL 01
 69 RCL 06
 70 RCL 05
 71 6
 72 ST+ Z
 73 +
 74 RCL IND Y
 75 RCL IND Y
 76 XEQ "K"
 77 RCL 06
 78 21
 79 +
 80 X<>Y
 81 RCL IND Y
 82 *
 83 RCL 05
 84 RCL 06              
 85 X=Y?
 86 GTO 00
 87 X<> Z
 88 GTO 01
 89 LBL 00
 90 RDN
 91 SIGN
 92 -
 93 LBL 01
 94 STO IND 04
 95 DSE 04
 96 DSE 05
 97 GTO 01
 98 15
 99 STO 05
100 DSE 06
101 GTO 01
102 276
103 STO 04
104 LBL 02
105 RCL 05
106 6
107 +
108 RCL IND X
109 XEQ "GG"
110 CHS
111 STO IND 04
112 DSE 04
113 DSE 05
114 GTO 02
115 37.27615
116 ENTER
117 CLX
118 XEQ "LS"
119 1/X
120 RCL 03             
121 LBL 10
122 STO 03
123 276
124 STO 04
125 21.006
126 STO 05
127 CLX
128 STO 06
129 LBL 03
130 RCL IND 05
131 RCL 03
132 XEQ "K"
133 RCL 05
134 15
135 +
136 X<>Y
137 RCL IND Y
138 *
139 RCL IND 04
140 *
141 ST+ 06
142 DSE 04
143 DSE 05
144 GTO 03
145 RCL 03
146 XEQ "GG"
147 ST+ 06
148 RCL 03             
149 SIGN
150 RCL 06
151 RTN
152 GTO 10
153 END

 
    ( 547 bytes / SIZE 277 )
 
 

      STACK        INPUTS      OUTPUTS
           X             x           f(x)
           L             /            x

 
Example:    Find f such that  f(x) = exp(-x) - x + §0+infinity x y f(y) dy
 
 

 01 LBL "GG"
 02 CHS
 03 E^X
 04 LASTX
 05 +
 06 RTN
 07 LBL "K"
 08 *
 09 RTN
 10 END

 
-If you want to evaluate f(0.7)

  0.7  XEQ "2FH0I+"  >>>>  0.496585304                                          ---Execution time = 17m32s---

-Likewise,

      0  R/S ( or XEQ 10 )  >>>>  1.000000000                                     ---Execution time = 15s---
      1  R/S ( or XEQ 10 )  >>>>  0.367879441

-The precision is again very good.
 

3°)  Fredholm Integral Equation 2nd kind over ]-infinity, +infinity[
 

  a) With Gauss-Hermite 10-point formula
 

 "2FHII" solves the integral equation:

   f(x) = g(x) + §-infinty+infinity K(x,y) f(y) dy    where a , b , g & K are known and f is the unknown function.

-We could also make the same change of argument ( y = Tan z ) and use "2FHAB"
-But "2FHII" employs Gauss-Hermite 10-point formula which leads to a linear system of 10 equations in 10 unknowns

-After solving this system, the integral equation is used again to evaluate f(x) for a given x
 

Data Registers:              R00 = det                R01-R02: unused

                                         R03 = x                   R11 to R30: abscissas & weights
                                         R04 to R06: temp    R31 to R140: coeff of the system

Flags: /
Subroutines:  "LS"  ( cf "Linear & Non-Linera Systems for the HP-41" )
                         "K"   which must take x & y in registers X & Y and return K(x,y) in X-register
                        "GG" ------------------  in X-register and return g(x) in X-register
 
 

 01 LBL "2FHII"
 02 STO 03
 03 130
 04 STO 04
 05 10
 06 STO 05
 07 STO 06
 08 .3429013272
 09 STO 11
 10 CHS
 11 STO 16
 12 1.03661083
 13 STO 12
 14 CHS
 15 STO 17
 16 1.756683649
 17 STO 13
 18 CHS
 19 STO 18
 20 2.532731674
 21 STO 14
 22 CHS
 23 STO 19
 24 3.436159119
 25 STO 15
 26 CHS
 27 STO 20
 28 .687081854
 29 STO 21
 30 STO 26
 31 .7032963231
 32 STO 22        
 33 STO 27
 34 .7414419319
 35 STO 23
 36 STO 28
 37 .8206661264
 38 STO 24
 39 STO 29
 40 1.025451691
 41 STO 25
 42 STO 30
 43 LBL 01
 44 RCL 06
 45 RCL 05
 46 10
 47 ST+ Z
 48 +
 49 RCL IND Y
 50 RCL IND Y
 51 XEQ "K"
 52 RCL 06        
 53 20
 54 +
 55 X<>Y
 56 RCL IND Y
 57 *
 58 RCL 05
 59 RCL 06
 60 X=Y?
 61 GTO 00
 62 X<> Z
 63 GTO 01
 64 LBL 00
 65 RDN
 66 SIGN
 67 -
 68 LBL 01
 69 STO IND 04
 70 DSE 04
 71 DSE 05
 72 GTO 01
 73 10
 74 STO 05        
 75 DSE 06
 76 GTO 01
 77 140
 78 STO 04
 79 LBL 02
 80 RCL 05
 81 10
 82 +
 83 RCL IND X
 84 XEQ "GG"
 85 CHS
 86 STO IND 04
 87 DSE 04
 88 DSE 05
 89 GTO 02
 90 31.1401
 91 ENTER
 92 CLX
 93 XEQ "LS"
 94 1/X
 95 RCL 03        
 96 LBL 10
 97 STO 03
 98 140
 99 STO 04
100 20.01
101 STO 05
102 CLX
103 STO 06
104 LBL 03
105 RCL IND 05
106 RCL 03
107 XEQ "K"
108 RCL 05
109 10
110 +
111 X<>Y
112 RCL IND Y
113 *
114 RCL IND 04
115 *
116 ST+ 06
117 DSE 04
118 DSE 05
119 GTO 03
120 RCL 03        
121 XEQ "GG"
122 ST+ 06
123 RCL 03
124 SIGN
125 RCL 06
126 RTN
127 GTO 10
128 END

 
    ( 313 bytes / SIZE 141 )
 
 

      STACK        INPUTS      OUTPUTS
           X             x           f(x)
           L             /            x

 
Example:    Find f such that   f(x) = exp(-x^2) - x sqrt(PI) + §-infinty+infinity 2 x.y2 f(y) dy
 
 

 01 LBL "GG"
 02 ENTER^
 03 X^2
 04 CHS
 05 E^X
 06 X<>Y
 07 PI
 08 SQRT
 09 *
 10 -
 11 RTN
 12 LBL "K"
 13 X<>Y
 14 X^2
 15 *
 16 ST+ X
 17 RTN
 18 END

 
-If you want to evaluate f(0.7)

  0.7  XEQ "2FHII"  >>>>  0.612626381                                          ---Execution time = 6m39s---

-Likewise,

      0  R/S ( or XEQ 10 )  >>>>  1.000000000                                          ---Execution time = 10s---
      1  R/S ( or XEQ 10 )  >>>>  0.367879422

-The precision is less good than before ( the exact solution is f(x) = exp(-x^2) )
  so,  f(0.7) = 0.612626394  &  f(1) = 0.367879441
 

  b) With Gauss-Hermite 15-point formula
 

 "2FHII+" solves the same integral equation but with Gauss-Hermite 15-point formula

   f(x) = g(x) + §-infinty+infinity K(x,y) f(y) dy    where a , b , g & K are known and f is the unknown function.
 

-After solving the 15x15 linear system, the integral equation is used again to evaluate f(x) for a given x
 

Data Registers:              R00 = det                R01-R02: unused

                                         R03 = x                   R07 to R36: abscissas & weights
                                         R04 to R06: temp    R37 to R276: coeff of the system

Flags: /
Subroutines:  "LS"  ( cf "Linear & Non-Linera Systems for the HP-41" )
                         "K"   which must take x & y in registers X & Y and return K(x,y) in X-register
                        "GG" ------------------  in X-register and return g(x) in X-register
 
 

 01 LBL "2FHII+"
 02 STO 03
 03 261
 04 STO 04
 05 15
 06 STO 05
 07 STO 06
 08 CLX
 09 STO 07
 10 .5650695833
 11 STO 08
 12 CHS
 13 STO 09
 14 1.136115585
 15 STO 10
 16 CHS
 17 STO 11
 18 1.719992575
 19 STO 12
 20 CHS
 21 STO 13
 22 2.325732486
 23 STO 14
 24 CHS
 25 STO 15
 26 2.967166928
 27 STO 16
 28 CHS
 29 STO 17
 30 3.669950373
 31 STO 18       
 32 CHS
 33 STO 19
 34 4.499990707
 35 STO 20
 36 CHS
 37 STO 21
 38 .5641003087
 39 STO 22
 40 .5670211534
 41 STO 23
 42 STO 24
 43 .5761933503
 44 STO 25
 45 STO 26
 46 .5930274498
 47 STO 27
 48 STO 28
 49 .6206626035
 50 STO 29
 51 STO 30
 52 .6661660051
 53 STO 31
 54 STO 32
 55 .748607366
 56 STO 33
 57 STO 34
 58 .9483689708
 59 STO 35
 60 STO 36
 61 LBL 01
 62 RCL 06
 63 RCL 05
 64 6
 65 ST+ Z
 66 +
 67 RCL IND Y
 68 RCL IND Y
 69 XEQ "K"
 70 RCL 06
 71 21
 72 +
 73 X<>Y
 74 RCL IND Y
 75 *
 76 RCL 05
 77 RCL 06
 78 X=Y?
 79 GTO 00
 80 X<> Z
 81 GTO 01
 82 LBL 00
 83 RDN
 84 SIGN
 85 -
 86 LBL 01
 87 STO IND 04
 88 DSE 04
 89 DSE 05
 90 GTO 01
 91 15
 92 STO 05       
 93 DSE 06
 94 GTO 01
 95 276
 96 STO 04
 97 LBL 02
 98 RCL 05
 99 6
100 +
101 RCL IND X
102 XEQ "GG"
103 CHS
104 STO IND 04
105 DSE 04
106 DSE 05
107 GTO 02
108 37.27615
109 ENTER
110 CLX
111 XEQ "LS"
112 1/X
113 RCL 03
114 LBL 10
115 STO 03
116 276
117 STO 04
118 21.006
119 STO 05
120 CLX
121 STO 06
122 LBL 03
123 RCL IND 05
124 RCL 03
125 XEQ "K"
126 RCL 05
127 15
128 +
129 X<>Y
130 RCL IND Y
131 *
132 RCL IND 04
133 *
134 ST+ 06
135 DSE 04
136 DSE 05
137 GTO 03
138 RCL 03
139 XEQ "GG"
140 ST+ 06
141 RCL 03
142 SIGN
143 RCL 06
144 RTN
145 GTO 10
146 END

 
    ( 389 bytes / SIZE 277 )
 
 

      STACK        INPUTS      OUTPUTS
           X             x           f(x)
           L             /            x

 
Example:    Find f such that   f(x) = exp(-x^2) - x sqrt(PI) + §-infinty+infinity 2 x.y2 f(y) dy
 
 

 01 LBL "GG"
 02 ENTER^
 03 X^2
 04 CHS
 05 E^X
 06 X<>Y
 07 PI
 08 SQRT
 09 *
 10 -
 11 RTN
 12 LBL "K"
 13 X<>Y
 14 X^2
 15 *
 16 ST+ X
 17 RTN
 18 END

 
-If you want to evaluate f(0.7)

  0.7  XEQ "2FHII+"  >>>>  0.612626392                                              ---Execution time = 17m11s---

-Likewise,

      0  R/S ( or XEQ 10 )  >>>>  1.000000000                                          ---Execution time = 15s---
      1  R/S ( or XEQ 10 )  >>>>  0.367879438

-The precision is better than with "2FHII"

-The exact results are:

       f(0.7) = 0.612626394  &  f(1) = 0.367879441   ( rounded to 9D )
 

4°) Volterra Integral Equation of the 1st kind
 

 "1VLTR" solves the integral equation:

     0 = g(x) + §ax K(x,y) f(y) dy    where a , g & K are known and f is the unknown function.

-The trapezoidal rule is used to approximate the integral and calculate  f(b) , assuming f(a) = g(a) = 0

-Let  xj = a + j.h    where   h = (b-a)/N  we have:

  0 = g0
  0 = g1 + §aa+h K(x,y) f(y) dy = g1 + (h/2) [ K(x1,x0) f(x0) + K(x1,x1) f(x1) ]  from which we can deduce  f(x1)

-We continue the same method until we find f(b)

-However, to increase the precision, the formula is applied with h , h/2 , h/4 , h/8 ( leading to 4 values, say A B C D )

-Then, a Romberg extrapolation to the limit is used with the 4 results obtained above and finally, f(b) is approximated by

         ( 4096 D - 1344 C + 84 B - A ) / 2835
 

Data Registers:              R00 = temp           ( Registers R01 to R03 are to be initialized before executing "1VLTR" )

                                      •  R01 = a                  R04 = h                   R16 = f(a)
                                      •  R02 = b                  R05 to R15: temp    R17 = f(a+h/8)       computed by the trapezoidal rule
                                      •  R03 = N                                                R18 = f(a+h/4)       with h = ( b - a ) / N
                                                                                                        .....................
Flags:   F02 is cleared by this program
Subroutines:    "K"   which must take x & y in registers X & Y and return K(x,y) in X-register
                         "GG" ------------------  in X-register and return g(x) in X-register
 
 
 

 01 LBL "1VLTR"
 02 CF 02
 03 GTO 14
 04 LBL "2VLTR"
 05 LBL 10
 06 SF 02
 07 LBL 14
 08 RCL 02        
 09 RCL 01
 10 -
 11 RCL 03
 12 /
 13 STO 04
 14 11.014
 15 STO 15
 16 LBL 00
 17 RCL 03
 18  E3
 19 /
 20 STO 05
 21 15
 22 STO 07
 23 LBL 01
 24 ISG 07
 25 CLX
 26 RCL 05
 27 INT
 28 RCL 04
 29 *
 30 RCL 01        
 31 +
 32 XEQ "GG"
 33 STO IND 07
 34 RCL 07
 35 STO 08
 36 LBL 02
 37 RCL 08
 38 16
 39 -
 40 X=0?
 41 GTO 04
 42 STO 09
 43 STO 10
 44 RCL 04
 45 *
 46 RCL 01
 47 +
 48 ENTER
 49 STO 06
 50 XEQ "K"
 51 STO 00
 52 LBL 03
 53 RCL 10        
 54 1
 55 ST- 08
 56 -
 57 STO 10
 58 RCL 04
 59 *
 60 RCL 01
 61 +
 62 RCL 06
 63 XEQ "K"
 64 RCL IND 08
 65 *
 66 RCL 04
 67 *
 68 1
 69 RCL 10
 70 X=0?
 71 ISG Y
 72 CLX
 73 RDN
 74 /
 75 ST+ IND 07
 76 RCL 10
 77 X#0?
 78 GTO 03
 79 CLX
 80 FS? 02
 81 SIGN
 82 RCL 00        
 83 RCL 04
 84 *
 85 2
 86 /
 87 -
 88 ST/ IND 07
 89 LBL 04
 90 ISG 05
 91 GTO 01
 92 RCL 03
 93 ST+ 03
 94 2
 95 ST/ 04
 96 RCL IND 07
 97 STO IND 15
 98 ISG 15
 99 GTO 00
100 16
101 ENTER
102 ST/ 03
103 ST* 04
104 RCL 03        
105 8
106 *
107 +
108  E3
109 /
110 +
111 RCL 14
112 4096
113 *
114 RCL 13
115 1344
116 *
117 -
118 RCL 12        
119 84
120 *
121 +
122 RCL 11
123 -
124 2835
125 /
126 FC? 02
127 GTO 14
128 RTN
129 GTO 10
130 LBL 14
131 END

 
     ( 203 bytes / SIZE var. )
 
 

      STACK        INPUTS      OUTPUTS
           Y             /    16.16+8.N
           X             /          f(b)

 
Example:          f(x) =  x exp(-x^2) - x + §0x 2 x.y f(y) dy     ->  Evaluate f(1)
 
 

 01 LBL "GG"
 02 ENTER^
 03 ENTER^
 04 X^2
 05 CHS
 06 E^X
 07 *
 08 X<>Y
 09 -
 10 RTN
 11 LBL "K"
 12 *
 13 ST+ X
 14 RTN
 15 END

 
    0  STO 01
    1  STO 02  and if you choose N = 2
    2  STO 03

    XEQ "1VLTR"  >>>>  f(1) = 0.367880022                                          ---Execution time = 4m20s---
                            X<>Y  16.032

 ( In R16 thru R32, you have f(0) f(1/16) f(1/8) ........... f(1)  computed by the trapezoidal rule but without Romberg extrapolation )

-The exact result is f(x) = exp(-x^2) thus,  f(1) = 0.367879441  ( rounded to 9D ) and the error = 58 E-8

-With N = 4 you will get   f(1) = 0.367879489
-And with N = 6 ,  f(1) = 0.367879448                 ( error = 7 E-9 )
 

5°) Volterra Integral Equation of the 2nd kind
 

 "2VLTR" solves the integral equation:

     f(x) = g(x) + §ax K(x,y) f(y) dy    where a , g & K are known and f is the unknown function.

-The trapezoidal rule is also used to approximate the integral.

-Let  xj = a + j.h    where   h = (b-a)/N  we have:

  f(x0) = g(x0)
  f(x1) = g(x1) + §aa+h K(x,y) f(y) dy = g1 + (h/2) [ K(x1,x0) f(x0) + K(x1,x1) f(x1) ]  from which we can deduce  f(x1)

-We continue the same method until we find f(b)

-However, to increase the precision, the formula is applied with h , h/2 , h/4 , h/8 ( leading to 4 values, say A B C D )

-Then, a Romberg extrapolation to the limit is used with the 4 results obtained above and finally, f(b) is approximated by

         ( 4096 D - 1344 C + 84 B - A ) / 2835
 

Data Registers:              R00 = temp             ( Registers R01 to R03 are to be initialized before executing "2VLTR" )

                                      •  R01 = a                  R04 = h                   R16 = f(a)
                                      •  R02 = b                  R05 to R15: temp    R17 = f(a+h/8)       computed by the trapezoidal rule
                                      •  R03 = N                                                R18 = f(a+h/4)       with h = ( b - a ) / N
                                                                                                        .....................
Flags:  F02 is set by this program
Subroutines:    "K"   which must take x & y in registers X & Y and return K(x,y) in X-register
                         "GG" ------------------  in X-register and return g(x) in X-register

>>>> The listing is given in the paragraph above
 
 

      STACK        INPUTS      OUTPUTS
           Y             /    16.16+8.N
           X             /          f(b)

 
Example:         f(x) = ( 1 + x ) exp(-x^2) - x + §0x 2 x.y f(y) dy     ->  Evaluate f(1)
 
 

 01 LBL "GG"
 02 ENTER^
 03 ENTER^
 04 X^2
 05 CHS
 06 E^X
 07 ST* Y
 08 +
 09 X<>Y
 10 -
 11 RTN
 12 LBL "K"
 13 *
 14 ST+ X
 15 RTN
 16 END

 
    0  STO 01
    1  STO 02  and if you choose N = 2
    2  STO 03

    XEQ "2VLTR"  >>>>  f(1) = 0.367880143                                          ---Execution time = 4m20s---
                            X<>Y  16.032

 ( In R16 thru R32, you have f(0) f(1/16) f(1/8) ........... f(1)  computed by the trapezoidal rule but without Romberg extrapolation )

-The exact result is f(x) = exp(-x^2) thus,  f(1) = 0.367879441  ( rounded to 9D ) and the error = 7 E-7

-With N = 4 you will get   f(1) = 0.367879444
-And with N = 8 ,  f(1) = 0.367879441                 ( error = -3 E-10 )

-But the execution time increases too...
 
 
 

Reference:

[1]  Press , Vetterling , Teukolsky , Flannery - "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4