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
 
 

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall