Poisson

# The Laplace and Poisson Equations for the HP-41

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
-The interval  [0,L'] is divided into N  parts

-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)
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
A program which computes   f2(x)   assuming x is in X-register upon entry
A program which computes    g1(y)  assuming y is in X-register upon entry
A program which computes    g2(y)  assuming y is in X-register upon entry

 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  XY 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"
E^X
ENTER^
1/X
-
2
/
RTN
LBL "U"
CLX
RTN
LBL "V"
CLX
RTN
LBL "W"
COS
1
E^X
ENTER^
1/X
-
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.
The numbers written in  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
-The interval  [0,L'] is divided into N parts

-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)
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.

A program which computes  f1(x)   assuming x is in X-register upon entry
A program which computes   f2(x)   assuming x is in X-register upon entry
A program which computes    g1(y)  assuming y is in X-register upon entry
A program which computes    g2(y)  assuming y is in X-register upon entry

-Line 173 is a three-byte  GTO 05

 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  XY 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  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.
The numbers written in  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.

A program which computes  f1(x)   assuming x is in X-register upon entry
A program which computes   f2(x)   assuming x is in X-register upon entry
A program which computes    g1(y)  assuming y is in X-register upon entry
A program which computes    g2(y)  assuming y is in X-register upon entry

-Line 195 is a three-byte  GTO 07

 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"    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  XY 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  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