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