The Diffusion Equation for the HP-41
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
A program which computes b(x,t) assuming x is in
X-register and t is in Y-register upon entry
A program which computes c(x,t) assuming x is in
X-register and t is in Y-register upon entry
A program which computes F(x) assuming x is
in X-register upon entry
A program which computes f(t)
assuming t is in X-register upon entry
A program which computes g(t) assuming
t is in X-register upon entry
-Line 102 is a three-byte GTO 02
-Line 111 is a three-byte GTO 01
| 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 103 RCL 15 104 INT 105 E3 106 / 107 20 108 + 109 RCL 00 110 RTN 111 GTO 01 112 END |
( 171 bytes / SIZE 21+M )
| STACK | INPUTS | OUTPUTS |
| Y | / | bbb.eee |
| X | / | t0+N.k |
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"
X^2 2 / RTN LBL "L" * CHS RTN LBL "M" SIGN CHS RTN LBL "N" X^2 ENTER^ SIGN + RTN LBL "O" CHS E^X RTN LBL "P" CHS E^X LASTX X^2 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
A program which computes b(x,t) assuming x
is in X-register and t is in Y-register upon entry
A program which computes c(x,t) assuming x
is in X-register and t is in Y-register upon entry
A program which computes F(x) assuming
x is in X-register upon entry
A program which computes f(t)
assuming t is in X-register upon entry
A program which computes g(t)
assuming t is in X-register upon entry
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
-Line 138 is a three-byte GTO 02
-Line 146 is a three-byte GTO 01
| 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 139 CLX 140 E3 141 / 142 20 143 + 144 RCL 00 145 RTN 146 GTO 01 147 END |
( 229 bytes / SIZE 4M+15 )
| STACK | INPUTS | OUTPUTS |
| Y | / | bbb.eee |
| X | / | t0+N.k |
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
¶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
A program which computes b(x,t) assuming x is in
X-register and t is in Y-register upon entry
A program which computes c(x,t) assuming x is in
X-register and t is in Y-register upon entry
A program which computes F(x) assuming x is
in X-register upon entry
A program which computes f(t)
assuming t is in X-register upon entry
A program which computes g(t) assuming
t is in X-register upon entry
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
-Line 158 is a three-byte GTO 02
-Line 166 is a three-byte GTO 01
| 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 159 RCL 17 160 E3 161 / 162 20 163 + 164 RCL 00 165 RTN 166 GTO 01 167 END |
( 260 bytes / SIZE 5M+15 )
| STACK | INPUTS | OUTPUTS |
| Y | / | bbb.eee |
| X | / | t0+N.k |
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