The Museum of HP Calculators

# The Diffusion Equation 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°) 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