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°) An Explicit Method
2°) An Implicit Method
3°) The Crank-Nicolson Scheme
-The following programs solve the partial differential equation ¶T/¶t = a(x,t) ¶2T/¶x2 + b(x,t) ¶T/¶x + c(x,t) T
or - using other notations: Tt = a(x,t) Txx + b(x,t) Tx + c(x,t) T where T , a , b , c are functions of x and t
with the initial values:
T(x,0) = F(x)
and the boundary conditions: T(0,t) = f(t)
a , b , c , F , f , g are known functions
T(L,t) = g(t)
-The interval [0,L] is divided into M parts.
x
L|-----------------------------------------------------
|
|
|
--|----------------------------------------------------- t
0
-The first routine uses an explicit method of order 1 in time t and
of order 2 in space x
-The second one uses an implicit method of the same order but it is
more stable.
-The third program employs an implicit - stable - method of order 2
in both space and time, thus producing more accurate results.
>>> The diffusion coefficient a(x,t) is assumed to be non-negative. Otherwise, the explicit method may be stable whereas the implicit methods are not!
-All these routines use the REGMOVE function from the X-Functions module.
1°) An Explicit Method
-The partial derivatives are approximated by ¶T/¶t
= Tt ~ [ T(x,t+k) - T(x,t) ]/k
where k = time stepsize
¶T/¶x
= Tx ~ [ T(x+h,t) - T(x-h,t) ]/(2h)
h = L/M
¶2T/¶x2
= Txx ~ [ T(x+h,t) - 2.T(x,t) + T(x+h,t) ]/h2
-The equation becomes Tm,n+1 =
Tm-1,n ( a.k/h2 - b.k/(2h) ) + Tm,n (
1 - 2.a.k/h2 + c.k ) + Tm+1,n ( a.k/h2
+ b.k ) where Tm,n =
T( m.h , n.k )
Data Registers: • R00 = t0 ( Registers R00 thru R10 are to be initialized before executing "DIF" )
• R01 = a name •
R04 = F name • R07
= L
• R10 = N = number of steps
• R02 = b name •
R05 = f name •
R08 = M
• R03 = c name •
R06 = g name • R09
= k
R11 = h = L/M , R12 to R18: temp R19: unused
>>>> When the program stops, R20 = T(0,t) R21 = T(h,t) R22 = T(2h,t) ................. R20+M = T(L,t) with t = t0 + N.k
Flags: /
Subroutines: A program which computes a(x,t)
assuming x is in X-register and t is in Y-register upon entry
--------------------------- b(x,t) ----------------------------------------------------------
--------------------------- c(x,t) ----------------------------------------------------------
--------------------------- F(x) ---------------------------
upon entry
--------------------------- f(t) --------- t ---------------------------
--------------------------- g(t) --------------------------------------
01 LBL "DIF"
02 RCL 07
03 RCL 08
04 STO 15
05 /
06 STO 11
07 LBL 00
08 RCL 11
09 RCL 15
10 *
11 XEQ IND 04
12 RCL 15
13 20
14 +
15 X<>Y
16 STO IND Y
17 DSE 15
18 GTO 00
19 CLX
20 XEQ IND 04
21 STO 20
22 LBL 01
23 RCL 10
24 STO 18
25 LBL 02
26 CLX
27 STO 14
28 20
29 STO 16
30 DSE X
31 RCL 08
32 +
33 E3
34 /
35 22
36 STO 17
37 DSE X
38 +
39 STO 15
40 LBL 03
41 RCL 11
42 ST+ 14
43 RCL 00
44 RCL 14
45 XEQ IND 01
46 RCL 11
47 X^2
48 /
49 STO 12
50 RCL 00
51 RCL 14
52 XEQ IND 02
53 RCL 11
54 ST+ X
55 /
56 STO 13
57 RCL 00
58 RCL 14
59 XEQ IND 03
60 RCL 12
61 ST+ X
62 -
63 RCL 09
64 ST* 12
65 ST* 13
66 *
67 RCL IND 15
68 ST* Y
69 +
70 RCL 13
71 ENTER^
72 CHS
73 RCL 12
74 ST+ Z
75 +
76 ST* IND 16
77 CLX
78 RCL IND 17
79 *
80 +
81 ST+ IND 16
82 1
83 ST+ 16
84 ST+ 17
85 ISG 15
86 GTO 03
87 19.02
88 RCL 08
89 E6
90 /
91 +
92 REGMOVE
93 RCL 09
94 ST+ 00
95 RCL 00
96 XEQ IND 05
97 STO 20
98 RCL 00
99 XEQ IND 06
100 STO IND 15
101 DSE 18
102 GTO 02
a three-byte GTO
103 RCL 15
104 INT
105 E3
106 /
107 20
108 +
109 RCL 00
110 RTN
111 GTO 01
a three-byte GTO
112 END
( 171 bytes / SIZE 21+M )
| STACK | INPUTS | OUTPUTS |
| Y | / | bbb.eee |
| X | / | t0+N.k |
where bbb.eee = control number of the solution bbb = 020 , eee = 20+M
Example: ¶T/¶t = (x2/2) ¶2T/¶x2 - t.x ¶T/¶x - T [ or Tt = (x2/2) Txx - t.x Tx - T ]
initial values:
T(x,0) = 1 + x2 ( t0
= 0 ) boundary conditions:
T(0,t) = exp(-t)
T(1,t) = exp(-t) + exp(-t2)
( L = 1 )
LBL "K"
LBL "L" LBL "M"
LBL "N" LBL "O"
LBL "P"
X^2
*
SIGN
X^2
CHS
CHS
2
CHS
CHS
ENTER^ E^X
E^X
/
RTN
RTN
SIGN
RTN
LASTX
RTN
+
X^2
RTN
CHS
E^X
+
RTN
t0 = 0 STO 00
"K" ASTO 01 "N"
ASTO 04 L = 1
STO 07
"L" ASTO 02
"O" ASTO 05
"M" ASTO 03 "P" ASTO
06
If we divide the interval [0,1] into M = 8
parts, 8
STO 08
If we choose the time stepsize
k = 1/32 32 1/X STO
09 and N = number of steps
= 2 STO 10
XEQ "DIF" >>>> t = 0.0625
( in 44 seconds )
X<>Y 20.028
and T(m/8 , 1/16) are in registers R20 to R28 ( m = 0 , 1 , ......... , 8 )
-Simply press R/S ( or XEQ 01 ) to continue with the same k and N (
or modify these parameters in R09 & R10 if needed )
-We find the following results:
| t\x | 0 | 1/8 | 2/8 | 3/8 | 4/8 | 5/8 | 6/8 | 7/8 | L=1 |
| 0 | 1 | 1.0156 | 1.0625 | 1.1406 | 1.2500 | 1.3906 | 1.5625 | 1.7656 | 2 |
| 1/16 | 0.9394 | 0.9541 | 1.0009 | 1.0788 | 1.1880 | 1.3283 | 1.4999 | 1.7022 | 1.9355 |
| 2/16 | 0.8825 | 0.8962 | 0.9425 | 1.0197 | 1.1278 | 1.2667 | 1.4364 | 1.6364 | 1.8670 |
| ... | ....... | ....... | ....... | ....... | ....... | ....... | ....... | ....... | ....... |
| 1 | 0.3679 | 0.3697 | 0.3861 | 0.4147 | 0.4550 | 0.5069 | 0.5718 | 0.6459 | 0.7358 |
| R20 | R21 | R22 | R23 | R24 | R25 | R26 | R27 | R28 |
-The values in black are the initial or boundary values.
-The numbers in red are computed by finite differencing.
-The exact solution is T(x,t) = exp(-t) + x2 exp(-t2)
Remarks
-If a , b , c are constants replace lines 41 to 73 by
RCL 01 RCL 02
RCL 03 *
RCL IND 15 ENTER^
RCL 11 RCL 11
RCL 09 R^
ST* Y
CHS
X^2
ST+ X ST* Z
ST+ X
+
R^
/
/
ST* T -
X<>Y
delete lines 26-27 and store the constants a , b , c
themselves into R01 R02 R03 ( instead of their
names )
The routine will run much faster!
-Numerical instability is a major problem with this method, if k is
not small enough.
-For instance, if we try to compute T(x,1) with k = 1/16, we get:
T(m/8,1) = 0.3656 ; 0.6414 ; -14.1373 ; 260.0787 ; -2055.3820 ; 7841.0783 ; -12672.4335 ( m = 1 , 2 , 3 , 4 , 5 , 6 , 7 )
-This oscillation is very far from the correct solution.
-I don't know the general rule but if b = c = 0 and if a is constant,
we must choose k <= h2/(2a)
-It usually requires very small k-values and execution time becomes
prohibitive.
-The 2 routines hereafter avoid this pitfall.
2°) An Implicit Method
-Here, the partial derivative with respect to t is approximated by
¶T/¶t
= Tt ~ [ T(x,t) - T(x,t-k) ]/k
with k = time stepsize
-The equation becomes Tm-1,n ( a.k/h2 - b.k/(2h) ) + Tm,n ( -1 - 2.a.k/h2 + c.k ) + Tm+1,n ( a.k/h2 + b.k ) = Tm,n-1 where Tm,n = T( m.h , n.k )
-So, the new values Tm-1,n ,
Tm,n , Tm+1,n
are the solutions of a tridiagonal system of M-1 equations
Data Registers: • R00 = t0 ( Registers R00 thru R10 are to be initialized before executing "DIF2" )
• R01 = a name •
R04 = F name • R07
= L
• R10 = N = number of steps
• R02 = b name •
R05 = f name •
R08 = M > 2
• R03 = c name •
R06 = g name • R09
= k
R11 = h = L/M , R12 to R16 & R18-R19: temp
R17: unused
R21+M to R4M+14: temp
>>>> When the program stops, R20 = T(0,t) R21 = T(h,t) R22 = T(2h,t) ................. R20+M = T(L,t) with t = t0 + N.k
Flag: F07 is cleared by this routine
Subroutines: A program which computes a(x,t)
assuming x is in X-register and t is in Y-register upon entry
--------------------------- b(x,t) ----------------------------------------------------------
--------------------------- c(x,t) ----------------------------------------------------------
--------------------------- F(x) ---------------------------
upon entry
--------------------------- f(t) --------- t ---------------------------
--------------------------- g(t) --------------------------------------
and "3DLS" tridiagonal systems ( cf "Linear and non-linear
systems for the HP-41" )
>>> delete lines 27-17 in the "3DLS" listing to avoid disturbing
R00
otherwise, add RCL 17 STO 00 X<>Y after
line 122 hereunder
and add RCL 00 STO 17 after line 116
01 LBL "DIF2"
02 RCL 07
03 RCL 08
04 STO 15
05 /
06 STO 11
07 LBL 00
08 RCL 11
09 RCL 15
10 *
11 XEQ IND 04
12 RCL 15
13 20
14 +
15 X<>Y
16 STO IND Y
17 DSE 15
18 GTO 00
19 CLX
20 XEQ IND 04
21 STO 20
22 LBL 01
23 RCL 10
24 STO 18
25 LBL 02
26 CLX
27 STO 14
28 RCL 08
29 3
30 *
31 16
32 +
33 STO 15
34 E3
35 /
36 21
37 +
38 STO 16
39 RCL 08
40 DSE X
41 E6
42 /
43 +
44 REGMOVE
45 RCL 09
46 ST+ 00
47 RCL 00
48 XEQ IND 05
49 STO 20
50 RCL 00
51 XEQ IND 06
52 STO 19
53 2 E-5
54 ST+ 16
55 SF 07
56 LBL 03
57 RCL 11
58 ST+ 14
59 RCL 00
60 RCL 14
61 XEQ IND 01
62 RCL 11
63 X^2
64 /
65 STO 12
66 RCL 00
67 RCL 14
68 XEQ IND 02
69 RCL 11
70 ST+ X
71 /
72 STO 13
73 RCL 12
74 -
75 RCL 09
76 ST* 12
77 ST* 13
78 *
79 FC?C 07
80 GTO 04
81 RCL 20
82 *
83 ST- IND 15
84 GTO 05
85 LBL 04
86 STO IND 16
87 ISG 16
88 LBL 05
89 RCL 00
90 RCL 14
91 XEQ IND 03
92 RCL 09
93 *
94 RCL 12
95 ST+ X
96 X<>Y
97 -
98 1
99 +
100 STO IND 16
101 RCL 12
102 RCL 13
103 +
104 CHS
105 ISG 16
106 FS? 30
107 GTO 06
108 STO IND 16
109 1
110 ST+ 15
111 ST- 16
112 GTO 03
113 LBL 06
114 RCL 19
115 *
116 ST- IND 15
117 RCL 15
118 E3
119 /
120 21
121 +
122 XEQ "3DLS"
123 INT
124 .021
125 +
126 RCL 08
127 DSE X
128 E6
129 /
130 +
131 REGMOVE
132 RCL 08
133 20
134 +
135 RCL 19
136 STO IND Y
137 DSE 18
138 GTO 02
a three-byte GTO
139 CLX
140 E3
141 /
142 20
143 +
144 RCL 00
145 RTN
146 GTO 01
a three-byte GTO
147 END
( 229 bytes / SIZE 4M+15 )
| STACK | INPUTS | OUTPUTS |
| Y | / | bbb.eee |
| X | / | t0+N.k |
where bbb.eee = control number of the solution bbb = 020 , eee = 20+M
Example: With the same partial differential equation, t0 = 0 STO 00 , ............ , k = 1/32 STO 09 & N = 2 STO 10
XEQ "DIF2" >>>> t = 0.0625
( in 68 seconds )
X<>Y 20.028
and T(m/8 , 1/16) are in registers R20 to R28 ( m = 0 , 1 , ......... , 8 )
-Simply press R/S ( or XEQ 01 ) to continue with the same k and N (
or modify these parameters in R09 & R10 if needed ), it gives:
| t\x | 0 | 1/8 | 2/8 | 3/8 | 4/8 | 5/8 | 6/8 | 7/8 | L=1 |
| 0 | 1 | 1.0156 | 1.0625 | 1.1406 | 1.2500 | 1.3906 | 1.5625 | 1.7656 | 2 |
| 1/16 | 0.9394 | 0.9558 | 1.0024 | 1.0801 | 1.1889 | 1.3287 | 1.4997 | 1.7019 | 1.9355 |
| 2/16 | 0.8825 | 0.8994 | 0.9455 | 1.0221 | 1.1294 | 1.2674 | 1.4363 | 1.6361 | 1.8670 |
| ... | ....... | ....... | ....... | ....... | ....... | ....... | ....... | ....... | ....... |
| 1 | 0.3679 | 0.3774 | 0.3955 | 0.4244 | 0.4646 | 0.5159 | 0.5784 | 0.6518 | 0.7358 |
| R20 | R21 | R22 | R23 | R24 | R25 | R26 | R27 | R28 |
-If we want to compute T(x,1) with k = 1/16, we get the following results:
T(m/8,1) = 0.3811 ; 0.4000 ; 0.4291 ; 0.4691 ; 0.5202 ; 0.5819 ; 0.6540 ( m = 1 , 2 , 3 , 4 , 5 , 6 , 7 )
-Of course, the accuracy is not so good as with k = 1/32
but instability is avoided.
3°) The Crank-Nicolson Scheme
-The approximation ¶T/¶t
= Tt ~ [ T(x,t+k) - T(x,t) ]/k
is of order 1
but this formula becomes second-order if it approximates the
derivative at t + k/2
-So, if we use the averages
¶T/¶x
= Tx ~ { [ T(x+h,t) - T(x-h,t) ]/(2h) +
[ T(x+h,t+k) - T(x-h,t+k) ]/(2h) } / 2
and ¶2T/¶x2
= Txx ~ { [ T(x+h,t) - 2.T(x,t) + T(x+h,t) ]/h2
+ [ T(x+h,t+k) - 2.T(x,t+k) + T(x+h,t+k) ]/h2 } /2
all these approximations are centered at t + k/2 and the formulas are of order 2 both in space and time.
-The diffusion equation becomes:
Tm-1,n+1 ( a/h2 - b/(2h) ) + Tm,n+1
( -2/k - 2a/h2 + c ) + Tm+1,n+1 ( a/h2
+ b/(2h) ) =
-Tm-1,n ( a/h2 - b/(2h) ) - Tm,n (
2/k - 2a/h2 + c ) - Tm+1,n ( a/h2 + b/(2h)
)
where Tm,n = T( m.h , n.k )
the functions a , b , c are to be evaluated at ( x , t + k/2 )
-Like "DIF2", "DIF3" must solve a tridiagonal linear system
of M-1 equations at each step.
Data Registers: • R00 = t0 ( Registers R00 thru R10 are to be initialized before executing "DIF3" )
• R01 = a name •
R04 = F name • R07
= L
• R10 = N = number of steps
• R02 = b name •
R05 = f name •
R08 = M > 2
• R03 = c name •
R06 = g name • R09
= k
R11 = h = L/M , R12 to R19: temp
R21+M to R5M+14: temp
>>>> When the program stops, R20 = T(0,t) R21 = T(h,t) R22 = T(2h,t) ................. R20+M = T(L,t) with t = t0 + N.k
Flag: F07 is cleared by this routine
Subroutines: A program which computes a(x,t)
assuming x is in X-register and t is in Y-register upon entry
--------------------------- b(x,t) ----------------------------------------------------------
--------------------------- c(x,t) ----------------------------------------------------------
--------------------------- F(x) ---------------------------
upon entry
--------------------------- f(t) --------- t ---------------------------
--------------------------- g(t) --------------------------------------
and "3DLS" tridiagonal systems ( cf "Linear and non-linear
systems for the HP-41" )
>>> delete lines 27-17 in the "3DLS" listing to avoid disturbing
R00
otherwise, add RCL 19 STO 00 X<>Y after
line 147 hereunder
and add RCL 00 STO 19 after line 139
01 LBL "DIF3"
02 RCL 07
03 RCL 08
04 STO 15
05 /
06 STO 11
07 LBL 00
08 RCL 11
09 RCL 15
10 *
11 XEQ IND 04
12 RCL 15
13 20
14 +
15 X<>Y
16 STO IND Y
17 DSE 15
18 GTO 00
19 CLX
20 XEQ IND 04
21 STO 20
22 LBL 01
23 RCL 10
24 STO 18
25 LBL 02
26 CLX
27 STO 14
28 RCL 08
29 4
30 *
31 16
32 +
33 STO 15
34 E3
35 /
36 21
37 +
38 RCL 08
39 +
40 2 E-5
41 +
42 STO 16
43 20
44 STO 17
45 RCL 09
46 ST+ 00
47 RCL 00
48 XEQ IND 05
49 STO 19
50 SF 07
51 RCL 09
52 2
53 /
54 ST- 00
55 LBL 03
56 CLX
57 STO IND 15
58 RCL 11
59 ST+ 14
60 RCL 00
61 RCL 14
62 XEQ IND 02
63 RCL 11
64 ST+ X
65 /
66 STO 13
67 RCL 00
68 RCL 14
69 XEQ IND 01
70 RCL 11
71 X^2
72 /
73 STO 12
74 RCL 13
75 -
76 FC?C 07
77 GTO 04
78 RCL 19
79 X<>Y
80 *
81 ST- IND 15
82 LASTX
83 GTO 05
84 LBL 04
85 STO IND 16
86 ISG 16
87 LBL 05
88 RCL IND 17
89 *
90 ST- IND 15
91 RCL 00
92 RCL 14
93 XEQ IND 03
94 RCL 12
95 ST+ 13
96 ST+ X
97 -
98 STO Y
99 RCL 09
100 1/X
101 ST+ X
102 ST+ Z
103 -
104 STO IND 16
105 X<>Y
106 ISG 17
107 CLX
108 RCL IND 17
109 *
110 ST- IND 15
111 ISG 17
112 CLX
113 RCL IND 17
114 RCL 13
115 *
116 ST- IND 15
117 ISG 16
118 FS? 30
119 GTO 06
120 LASTX
121 STO IND 16
122 1
123 ST+ 15
124 ST- 16
125 ST- 17
126 GTO 03
127 LBL 06
128 RCL 09
129 2
130 /
131 ST+ 00
132 RCL 00
133 XEQ IND 06
134 STO IND 17
135 RCL 13
136 *
137 ST- IND 15
138 RCL 19
139 STO 20
140 RCL 15
141 E3
142 /
143 21
144 +
145 RCL 08
146 +
147 XEQ "3DLS"
148 INT
149 .021
150 +
151 RCL 08
152 DSE X
153 E6
154 /
155 +
156 REGMOVE
157 DSE 18
158 GTO 02
a three-byte GTO
159 RCL 17
160 E3
161 /
162 20
163 +
164 RCL 00
165 RTN
166 GTO 01
a three-byte GTO
167 END
( 260 bytes / SIZE 5M+15 )
| STACK | INPUTS | OUTPUTS |
| Y | / | bbb.eee |
| X | / | t0+N.k |
where bbb.eee = control number of the solution bbb = 020 , eee = 20+M
Example: With the same equation, t0 = 0 STO 00 , ............ , k = 1/16 STO 09 & N = 1 STO 10
XEQ "DIF3" >>>> t = 0.0625
( in 41 seconds )
X<>Y 20.028
and T(m/8 , 1/16) are in registers R20 to R28 ( m = 0 , 1 , ......... , 8 )
-Simply press R/S ( or XEQ 01 ) to continue with the same k and N (
or modify these parameters in R09 & R10 if needed ), it yields:
| t\x | 0 | 1/8 | 2/8 | 3/8 | 4/8 | 5/8 | 6/8 | 7/8 | L=1 |
| 0 | 1 | 1.0156 | 1.0625 | 1.1406 | 1.2500 | 1.3906 | 1.5625 | 1.7656 | 2 |
| 1/16 | 0.9394 | 0.9550 | 1.0017 | 1.0795 | 1.1884 | 1.3285 | 1.4997 | 1.7020 | 1.9355 |
| 2/16 | 0.8825 | 0.8978 | 0.9440 | 1.0209 | 1.1286 | 1.2670 | 1.4362 | 1.6362 | 1.8670 |
| ... | ....... | ....... | ....... | ....... | ....... | ....... | ....... | ....... | ....... |
| 1 | 0.3679 | 0.3735 | 0.3908 | 0.4195 | 0.4597 | 0.5114 | 0.5747 | 0.6494 | 0.7358 |
| R20 | R21 | R22 | R23 | R24 | R25 | R26 | R27 | R28 |
-The maximum error is smaller than 2 E-4
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