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