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) + §_a^b 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⁺¹ ( -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) + §_a^b 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

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⁺¹ ( -x - 2.y ) exp(-x.y) f(y) dy

-Load these 2 subroutines in main memory

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) + §₀^+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

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 + §₀^+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) + §₀^+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

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 + §₀^+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

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.y² 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

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.y² 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) + §_a^x 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 x_j = a + j.h where h = (b-a)/N we have:

0 = g₀
0 = g₁ + §_a^a+h K(x,y) f(y) dy = g₁ + (h/2) [ K(x₁,x₀) f(x₀) + K(x₁,x₁) f(x₁) ] from which we can deduce f(x₁)

-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 + §₀^x 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) + §_a^x 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 x_j = a + j.h where h = (b-a)/N we have:

f(x₀) = g(x₀)
f(x₁) = g(x₁) + §_a^a+h K(x,y) f(y) dy = g₁ + (h/2) [ K(x₁,x₀) f(x₀) + K(x₁,x₁) f(x₁) ] from which we can deduce f(x₁)

-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 + §₀^x 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

hp41programs

Integral Equations for the HP-41