The Museum of HP Calculators
The Museum of HP Calculators
This program is Copyright © 2007 by Jean-Marc Baillard and is used here by permission.
This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.
1°) Laplace Equation ¶2u/¶x2
+ ¶2u/¶y2
= uxx + uyy = 0
2°) Poisson Equation ¶2u/¶x2
+ ¶2u/¶y2
= uxx + uyy = F(x,y)
a) Program#1
b) Program#2 ( X-Functions
Module required )
-The 3 programs hereafter use finite differencing to solve these partial
differential equations
with boundary conditions defined on a rectangle [0,L]x[0,L']
-They employ a method of order 2.
y
L'|----------------------|-
|
|
|
|
|
|
--|----------------------|--------x
0
L
1°) Laplace Equation
-This program solves the partial differential equation: ¶2u/¶x2 + ¶2u/¶y2 = uxx + uyy = 0
with the boundary conditions: u(x,0)
= f1(x) u(0,y) = g1(y)
0 <= x <= L
u(x,L') = f2(x) u(L,y) = g2(y)
0 <= y <= L'
-The interval [0,L] is divided into M parts
------------ [0,L'] -------------- N ----
-The partial derivatives are approximated by ¶2u/¶x2
= uxx ~ ( ui-1,j - 2.ui,j + ui+1,j
)/h2 where h = L/M
, ui,j = u(i.h,j.k)
and ¶2u/¶y2
= uyy ~ ( ui,j-1 - 2.ui,j + ui,j+1
)/k2 where k = L'/N
-It yields: ui-1,j + ui+1,j + h2/k2 ( ui,j-1 + ui,j+1 ) = 2.( 1 + h2/k2 ).ui,j
-So, we have to solve a linear system of (M-1).(N-1) equations in (M-1).(N-1)
unknowns.
-Fortunately, this is a sparse system and we can use an overrelaxation
method.
-The overrelaxation parameter is 1.2 ( line 02 ) but it's not
always the best choice:
change line 02 as you like or delete lines 02-03 and initialize
register R09 before executing this routine.
-The overrelaxation parameter is usually between 1 and 2.
-The iterations stop when the last correction is rounded to 0 ( line
156-157 )
-So, the accuracy is controlled by the display format.
-However, for instance, a 4-place accuracy in the solution of the linear
system
doesn't necessarily mean a 4-place accuracy in u(x,y)
Data Registers: ( Registers R01 thru R08 are to be initialized before executing "LAP" )
• R01 = f1 name •
R05 = L
R09 = 1.2
R00 & R11 to R16: temp
• R02 = f2 name •
R06 = M R10
= h2/k2
R17 to Rff = ui,j in column order ff = 16+(M+1).(N+1)
• R03 = g1 name •
R07 = L'
R11 = h = L/M
• R04 = g2 name •
R08 = N
R12 = k = L'/N
Flags: /
Subroutines: A program which computes
f1(x) assuming x is in X-register upon entry
--------------------------- f2(x)
-------------------------------------
--------------------------- g1(y) assuming
y is in X-register upon entry
--------------------------- g2(y) -------------------------------------
01 LBL "LAP"
02 1.2
03 STO 09
04 RCL 05
05 RCL 06
06 /
07 STO 11
08 LASTX
09 STO 14
10 RCL 08
11 1
12 ST+ Z
13 +
14 *
15 16
16 +
17 STO 13
18 STO 15
19 LBL 02
20 RCL 11
21 RCL 14
22 *
23 XEQ IND 02
24 STO IND 13
25 DSE 13
26 DSE 14
27 GTO 02
28 RCL 06
29 E-5
30 ST* Y
31 +
32 ST+ 13
33 ST+ 15
34 RCL 07
35 RCL 08
36 STO 10
37 STO 14
38 /
39 STO 12
40 LBL 03
41 RCL 12
42 RCL 14
43 *
44 XEQ IND 03
45 STO IND 13
46 DSE 13
47 DSE 14
48 GTO 03
49 LBL 04
50 RCL 10
51 RCL 12
52 *
53 XEQ IND 04
54 STO IND 15
55 DSE 15
56 DSE 10
57 GTO 04
58 RCL 15
59 INT
60 STO 13
61 RCL 06
62 STO 14
63 LBL 01
64 RCL 11
65 RCL 14
66 *
67 XEQ IND 01
68 STO IND 13
69 DSE 13
70 DSE 14
71 GTO 01
72 CLX
73 XEQ IND 01
74 STO 17
75 RCL 11
76 RCL 12
77 /
78 X^2
79 STO 10
80 LBL 05
81 18
82 STO 12
83 RCL 06
84 STO 16
85 +
86 STO 13
87 1
88 +
89 STO 14
90 LASTX
91 +
92 STO 15
93 RCL 06
94 RCL 08
95 LASTX
96 ST+ Z
97 +
98 *
99 15
100 +
101 E3
102 /
103 +
104 ST+ 16
105 CLX
106 STO 00
107 LBL 06
108 RCL 06
109 STO 11
110 DSE 11
111 LBL 07
112 RCL IND 12
113 RCL IND 16
114 +
115 RCL 10
116 *
117 RCL IND 13
118 RCL IND 15
119 +
120 +
121 RCL 10
122 ENTER^
123 SIGN
124 +
125 ST+ X
126 /
127 RCL IND 14
128 -
129 RCL 09
130 *
131 ST+ IND 14
132 ABS
133 RCL 00
134 X<Y?
135 X<>Y
136 STO 00
137 SIGN
138 ST+ 12
139 ST+ 13
140 ST+ 14
141 ST+ 15
142 ISG 16
143 X<0?
144 GTO 08
145 DSE 11
146 GTO 07
147 2
148 ST+ 12
149 ST+ 13
150 ST+ 14
151 ST+ 15
152 ST+ 16
153 GTO 06
154 LBL 08
155 LASTX
156 RND
157 X#0?
158 GTO 05
159 RCL 16
160 INT
161 E3
162 /
163 17
164 +
165 RCL 06
166 1
167 +
168 E5
169 /
170 +
171 END
( 243 bytes / SIZE 17+(M+1)(N+1) )
| STACK | INPUT | OUTPUT |
| X | / | bb.eee,ii |
where bb.eee,ii is the control number of the array: bb = 17 , eee = 16+(M+1)(N+1) , ii = N+1
Example:
¶2u/¶x2 + ¶2u/¶y2 = uxx + uyy = 0 0 <= x <= 1 , 0 <= y <= PI/2 ( so L = 1 , L' = PI/2 )
boundary conditions:
u(x,0) = sinh x
u(0,y) = 0
u(x,pi/2) = 0
u(1,y) = (sinh 1) cos y
LBL "T" LBL "U"
LBL "V" LBL "W"
E^X
CLX
CLX COS
ENTER^ RTN
RTN 1
1/X
E^X
-
ENTER^
2
1/X
/
-
RTN
2
/
*
RTN
"T" ASTO 01
"U" ASTO 02
"V" ASTO 03
here, we could delete LBL "V" and store "U" in R03
"W" ASTO 04
L = 1 STO 05 L' = PI/2 STO 07 XEQ RAD
-If ui,j = 0 is your initial guess, clear registers
R17 to R51 17.051 XEQ "CLRGX" ( or SIZE
017 , SIZE 052 - or greater )
-If we choose FIX 3 , M = 4 and
N = 6
FIX 3
4 STO 06
6 STO 08
XEQ "LAP" >>>> 17.05105
( in 3mn45s ) and we have the following results in registers
R17 thru R51 ( rounded to 4 decimals )
R17 R22
R27 R32
R37 R42
R47
R18 R23
R28
R33
R38
R43
R48
R19 R24
R29
R34
R39
R44
R49
R20 R25
R30
R35
R40
R45
R50
R21 R26
R31 R36
R41 R46
R51
| x \ y | 0 | PI/12 | 2PI/12 | 3PI/12 | 4PI/12 | 5PI/12 | PI/2 |
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 1/4 | 0.2526 | 0.2441 | 0.2189 | 0.1788 | 0.1264 | 0.0655 | 0 |
| 2/4 | 0.5211 | 0.5036 | 0.4516 | 0.3688 | 0.2608 | 0.1350 | 0 |
| 3/4 | 0.8223 | 0.7946 | 0.7125 | 0.5818 | 0.4114 | 0.2130 | 0 |
| 1 | 1.1752 | 1.1352 | 1.0178 | 0.8310 | 0.5876 | 0.3042 | 0 |
-The numbers written in black are the boundary-values.
----------------------- red are computed by
the finite difference method
and they approximate the solution of a linear system of
15 equations in 15 unknowns.
-For instance, u(3/4,PI/12) ~ 0.7946 whereas the correct value is 0.794297
-The exact solution is u(x,y) = sinh x cos y
Note:
-We can get a better accuracy with larger M- and N-values and if we
execute "LAP" in FIX 6 or greater
-With M = 8 , N = 12 it gives u(3/4,PI/12) ~ 0.794380
= R41
2°) Poisson Equation
a) Program#1
-The following program solves the partial differential equation: ¶2u/¶x2 + ¶2u/¶y2 = uxx + uyy = F(x,y) where F is a source term
with the boundary conditions: u(x,0)
= f1(x) u(0,y) = g1(y)
0 <= x <= L
u(x,L') = f2(x) u(L,y) = g2(y)
0 <= y <= L'
-The interval [0,L] is divided into M parts
------------ [0,L'] -------------- N ----
-The partial derivatives are approximated by ¶2u/¶x2
= uxx ~ ( ui-1,j - 2.ui,j + ui+1,j
)/h2 where h = L/M
, ui,j = u(i.h,j.k)
and ¶2u/¶y2
= uyy ~ ( ui,j-1 - 2.ui,j + ui,j+1
)/k2 where k = L'/N
-It yields: -h2.Fi,j + ui-1,j + ui+1,j + h2/k2 ( ui,j-1 + ui,j+1 ) = 2.( 1 + h2/k2 ).ui,j where Fi,j = F(i.h,j.k)
-We use the same overrelaxation method.
-The overrelaxation parameter is 1.2 ( line 02 ) but you can
change line 02 as you like
or delete lines 02-03 and initialize register R09 before executing
this routine.
-The best overrelaxation parameter is usually between 1 and 2.
-The iterations stop when the last correction is rounded to 0 ( line
171-172 )
-Thus, the accuracy is controlled by the display format.
-However, an accurate solution to the linear system is not always an
accurate solution to u(x,y)
Data Registers: • R00 = F name ( Registers R00 thru R08 are to be initialized before executing "POI" )
• R01 = f1 name •
R05 = L
R09 = 1.2
R13 to R21: temp
• R02 = f2 name •
R06 = M R10
= h2/k2
R22 to Rff = ui,j in column order ff = 21+(M+1).(N+1)
• R03 = g1 name •
R07 = L'
R11 = h = L/M
• R04 = g2 name •
R08 = N
R12 = k = L'/N
Flags: /
Subroutines: A program which computes F(x,y)
assuming x is in X-register and y is in Y-register upon
entry.
--------------------------- f1(x)
-------------------------- upon entry
--------------------------- f2(x)
-------------------------------------
--------------------------- g1(y) assuming
y is in X-register upon entry
--------------------------- g2(y) -------------------------------------
01 LBL "POI"
02 1.2
03 STO 09
04 RCL 05
05 RCL 06
06 /
07 STO 11
08 LASTX
09 STO 14
10 RCL 08
11 1
12 ST+ Z
13 +
14 *
15 21
16 +
17 STO 13
18 STO 15
19 LBL 02
20 RCL 11
21 RCL 14
22 *
23 XEQ IND 02
24 STO IND 13
25 DSE 13
26 DSE 14
27 GTO 02
28 RCL 06
29 E-5
30 ST* Y
31 +
32 ST+ 13
33 ST+ 15
34 RCL 07
35 RCL 08
36 STO 10
37 STO 14
38 /
39 STO 12
40 LBL 03
41 RCL 12
42 RCL 14
43 *
44 XEQ IND 03
45 STO IND 13
46 DSE 13
47 DSE 14
48 GTO 03
49 LBL 04
50 RCL 10
51 RCL 12
52 *
53 XEQ IND 04
54 STO IND 15
55 DSE 15
56 DSE 10
57 GTO 04
58 RCL 15
59 INT
60 STO 13
61 RCL 06
62 STO 14
63 LBL 01
64 RCL 11
65 RCL 14
66 *
67 XEQ IND 01
68 STO IND 13
69 DSE 13
70 DSE 14
71 GTO 01
72 CLX
73 XEQ IND 01
74 STO 22
75 RCL 11
76 RCL 12
77 /
78 X^2
79 STO 10
80 LBL 05
81 23
82 STO 16
83 RCL 06
84 STO 20
85 +
86 STO 17
87 1
88 +
89 STO 18
90 LASTX
91 +
92 STO 19
93 RCL 06
94 RCL 08
95 LASTX
96 ST+ Z
97 +
98 *
99 20
100 +
101 E3
102 /
103 +
104 ST+ 20
105 CLX
106 STO 14
107 STO 15
108 LBL 06
109 CLX
110 STO 13
111 RCL 12
112 ST+ 14
113 RCL 06
114 STO 21
115 DSE 21
116 LBL 07
117 RCL 11
118 ST+ 13
119 RCL 14
120 RCL 13
121 XEQ IND 00
122 RCL 11
123 X^2
124 *
125 RCL IND 16
126 RCL IND 20
127 +
128 RCL 10
129 *
130 X<>Y
131 -
132 RCL IND 17
133 RCL IND 19
134 +
135 +
136 RCL 10
137 ENTER^
138 SIGN
139 +
140 ST+ X
141 /
142 RCL IND 18
143 -
144 RCL 09
145 *
146 ST+ IND 18
147 ABS
148 RCL 15
149 X<Y?
150 X<>Y
151 STO 15
152 SIGN
153 ST+ 16
154 ST+ 17
155 ST+ 18
156 ST+ 19
157 ISG 20
158 X<0?
159 GTO 08
160 DSE 21
161 GTO 07
162 2
163 ST+ 16
164 ST+ 17
165 ST+ 18
166 ST+ 19
167 ST+ 20
168 GTO 06
169 LBL 08
170 LASTX
171 RND
172 X#0?
173 GTO 05
a three-byte GTO ( or replace this line by GTO 15 and line 80 by
LBL 15 )
174 RCL 20
175 INT
176 E3
177 /
178 22
179 +
180 RCL 06
181 1
182 +
183 E5
184 /
185 +
186 END
( 267 bytes / SIZE 22+(M+1)(N+1) )
| STACK | INPUT | OUTPUT |
| X | / | bb.eee,ii |
where bb.eee,ii is the control number of the array: bb = 22 , eee = 21+(M+1)(N+1) , ii = N+1
Example:
¶2u/¶x2 + ¶2u/¶y2 = uxx + uyy = 2.exp(-x-y) ; 0 <= x <= 1 , 0 <= y <= 1 ( L = L' = 1 )
boundary conditions:
u(x,0) = exp(-x)
u(0,y) = exp(-y)
u(x,1) = exp(-1-x)
u(1,y) = exp(-1-y)
LBL "T" LBL "U"
LBL "V" LBL "W"
LBL "Z"
CHS
1
CHS 1
+
E^X
+
E^X +
CHS
RTN
CHS
RTN CHS
E^X
E^X
E^X
ST+ X
RTN
RTN
RTN
"Z" ASTO 00
"T" ASTO 01
"U" ASTO 02
"V" ASTO 03
( "V" & "W" are actually unuseful here, they may be replaced by "T"
& "U" )
"W" ASTO 04
L = 1 STO 05 L' = 1 STO 07
-If ui,j = 0 is your initial guess, clear registers
R22 to R51 22.051 XEQ "CLRGX" ( or SIZE
022 , SIZE 052 - or greater )
-If we choose FIX 3 , M = 4 and
N = 5
FIX 3
4 STO 06
5 STO 08
XEQ "POI" >>>> 22.05105
( in 4mn06s ) and we have the following results in registers
R22 thru R51 ( rounded to 4 decimals )
R22 R27
R32 R37
R42 R47
R23 R28
R33
R38
R43
R48
R24 R29
R34
R39
R44
R49
R25 R30
R35
R40
R45
R50
R26 R31
R36 R41
R46 R51
| x \ y | 0 | 1/5 | 2/5 | 3/5 | 4/5 | 1 |
| 0 | 1 | 0.8187 | 0.6703 | 0.5488 | 0.4493 | 0.3679 |
| 1/4 | 0.7788 | 0.6377 | 0.5222 | 0.4276 | 0.3500 | 0.2865 |
| 2/4 | 0.6065 | 0.4967 | 0.4067 | 0.3330 | 0.2727 | 0.2231 |
| 3/4 | 0.4724 | 0.3868 | 0.3168 | 0.2594 | 0.2123 | 0.1738 |
| 1 | 0.3679 | 0.3012 | 0.2466 | 0.2019 | 0.1653 | 0.1353 |
-The numbers written in black are the boundary-values.
----------------------- red are computed by
the finite difference method
and they are the approximate solution of a linear system
of 12 equations in 12 unknowns.
-For instance, u(3/4,2/5) ~ 0.3168 whereas
the correct value is 0.316637
-The results are more accurate than what we could expect!
-The exact solution is u(x,y) = exp(-x-y)
Notes:
-We can get a better accuracy with larger M- and N-values and if we
execute "POI" in FIX 6
-For example, with M = 8 and N = 10 we find that register R64 = u(3/4,2/5)
~ 0.316677
-With M = 8 and N = 10, we solve a system of 63 equations.
-So, a good emulator - like Warren Furlow's V41 in turbo mode - is
quite useful if you want to execute this routine for ( relatively ) large
M- & N-values.
b) Program#2 ( X-Functions
Module required )
-"POI" re-computes F(x,y) at each iteration.
-The following routine uses a Data-file to store these values - thus
saving execution time.
Data Registers: • R00 = F name ( Registers R00 thru R08 are to be initialized before executing "POI2" )
• R01 = f1 name •
R05 = L
R09 = 1.2
• R02 = f2 name •
R06 = M R10
= h2/k2
h = L/M , k = L'/N
• R03 = g1 name •
R07 = L'
R11 to R17: temp
• R04 = g2 name •
R08 = N
R18 to Rff = ui,j in column order , ff = 17+(M+1).(N+1)
Flags: /
Subroutines: A program which computes F(x,y)
assuming x is in X-register and y is in Y-register upon
entry.
--------------------------- f1(x)
-------------------------- upon entry
--------------------------- f2(x)
-------------------------------------
--------------------------- g1(y) assuming
y is in X-register upon entry
--------------------------- g2(y) -------------------------------------
01 LBL "POI2"
02 1.2
03 STO 09
04 RCL 05
05 RCL 06
06 /
07 STO 11
08 LASTX
09 STO 14
10 RCL 08
11 1
12 ST+ Z
13 +
14 *
15 17
16 +
17 STO 13
18 STO 15
19 LBL 02
20 RCL 11
21 RCL 14
22 *
23 XEQ IND 02
24 STO IND 13
25 DSE 13
26 DSE 14
27 GTO 02
28 RCL 06
29 E-5
30 ST* Y
31 +
32 ST+ 13
33 ST+ 15
34 RCL 07
35 RCL 08
36 STO 10
37 STO 14
38 /
39 STO 12
40 LBL 03
41 RCL 12
42 RCL 14
43 *
44 XEQ IND 03
45 STO IND 13
46 DSE 13
47 DSE 14
48 GTO 03
49 LBL 04
50 RCL 10
51 RCL 12
52 *
53 XEQ IND 04
54 STO IND 15
55 DSE 15
56 DSE 10
57 GTO 04
58 RCL 15
59 INT
60 STO 13
61 RCL 06
62 STO 14
63 LBL 01
64 RCL 11
65 RCL 14
66 *
67 XEQ IND 01
68 STO IND 13
69 DSE 13
70 DSE 14
71 GTO 01
72 CLX
73 XEQ IND 01
74 STO 18
75 RCL 06
76 RCL 08
77 1
78 ST- Z
79 -
80 STO 10
81 *
82 "TEMP"
the size of this temporary data file is (M-1).(N-1)
83 CRFLD
84 CLX
85 STO 14
86 LBL 05
87 CLX
88 STO 13
89 RCL 12
90 ST+ 14
91 RCL 06
92 STO 15
93 DSE 15
94 LBL 06
95 RCL 11
96 ST+ 13
97 RCL 14
98 RCL 13
99 XEQ IND 00
100 RCL 11
101 X^2
102 *
103 SAVEX
104 DSE 15
105 GTO 06
106 DSE 10
107 GTO 05
108 RCL 11
109 RCL 12
110 /
111 X^2
112 STO 10
113 LBL 07
114 CLX
115 SEEKPT
116 19
117 STO 13
118 RCL 06
119 STO 17
120 +
121 STO 14
122 1
123 +
124 STO 15
125 LASTX
126 +
127 STO 16
128 RCL 06
129 RCL 08
130 LASTX
131 ST+ Z
132 +
133 *
134 16
135 +
136 E3
137 /
138 +
139 ST+ 17
140 CLX
141 STO 11
142 LBL 08
143 RCL 06
144 STO 12
145 DSE 12
146 LBL 09
147 RCL IND 13
148 RCL IND 17
149 +
150 RCL 10
151 *
152 RCL IND 14
153 RCL IND 16
154 +
155 +
156 GETX
157 -
158 RCL 10
159 ENTER^
160 SIGN
161 +
162 ST+ X
163 /
164 RCL IND 15
165 -
166 RCL 09
167 *
168 ST+ IND 15
169 ABS
170 RCL 11
171 X<Y?
172 X<>Y
173 STO 11
174 SIGN
175 ST+ 13
176 ST+ 14
177 ST+ 15
178 ST+ 16
179 ISG 17
180 X<0?
181 GTO 10
182 DSE 12
183 GTO 09
184 2
185 ST+ 13
186 ST+ 14
187 ST+ 15
188 ST+ 16
189 ST+ 17
190 GTO 08
191 LBL 10
192 LASTX
193 RND
194 X#0?
195 GTO 07
a three-byte GTO
196 RCL 17
197 INT
198 E3
199 /
200 18
201 +
202 RCL 06
203 1
204 +
205 E5
206 /
207 +
208 "TEMP"
209 PURFL
210 CLA
211 END
( 308 bytes / SIZE 18+(M+1)(N+1) )
| STACK | INPUT | OUTPUT |
| X | / | bb.eee,ii |
where bb.eee,ii is the control number of the array: bb = 18 , eee = 17+(M+1)(N+1) , ii = N+1
-The same example is solved in 3mn18s instead of 4mn06s
-The results are in R18 thru R47 instead of R22 to R51 ( the
output is 18.04705 )
References:
[1] Press , Vetterling , Teukolsky , Flannery , "Numerical Recipes
in C++" - Cambridge University Press - ISBN 0-521-75033-4
[2] F. Scheid - "Numerical Analysis" - Mc Graw Hill
ISBN 07-055197-9
Go back to the HP-41 software library
Go back to the general software library
Go
back to the main exhibit hall