 The Museum of HP Calculators

# The Laplace and Poisson Equations for the HP-41

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.

## Overview

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:

  Press , Vetterling , Teukolsky , Flannery , "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4
  F. Scheid - "Numerical Analysis" - Mc Graw Hill   ISBN 07-055197-9