Diffusion

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:

  Press , Vetterling , Teukolsky , Flannery , "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4
  F. Scheid - "Numerical Analysis" - Mc Graw Hill   ISBN 07-055197-9