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:

[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