Partial Differential Equations Module


Overview
 

-The following routines employs finite differences to solve a few partial differential equations:

  Generalized Diffusion equations: 1D & 2D, linear & non-linear
  2-Dimensional generalized Poisson equation
  Generalized wave equations: 1D & 2D

+ Euler-Lagrange Equation that appear in variational problems.

-In most cases, the partial derivatives are approximated by formulas of order 2:

    f/x  ~   [ f(x+h) - f(x-h) ] / 2h
   2f/x2 ~ [ f(x+h) - 2 f(x) + f(x+h) ] / h2

  2f/xy ~ [ f(x+h,y+k) + f(x-h,y-k) - f(x+h,y-k) - f(x-h,y+k) ] / (4.h.k)

-However, for non-liuear "diffusion equatrions"   f/x  ~   [ f(x+h) - f(x) ] / h  which is only 1st order
 
 
 
XROM  Function  Desciption
 19,00
 19,01
 19,02
 19,03
 19,04
 19,05
 19,06
 19,07
 19,08
 19,09
 19,10
 -PARDIFFEQ
  "3DLS"
  "DIFF"
  "2DIFF"
  "NLDIFF"
  "NL2DIFF"
  "2DPOIS"
  "WAVE"
  "2DWAVE"
  "SHOOT"
  "EULAG"
 Section Header
 Tridiagonal Linear System
 1 Dimensional Diffusion Equation ( Crank-Nicolson )
 2 Dimensional Diffusion Equation ( C-N )
 Non-Linear Diffusion Equation ( explicit )
 Non-Linear 2D Diffusion Equation ( explicit )
 2D Poisson Equation
 Wave Equation ( explicit )
 2D Wave Equation ( explicit )
 y"=T(x,y,y') knowing y(a) & y(b) ( garden-hose )
 Variational Problems: Euler-Lagrange Equation.

 

Tridiagonal Linear Systems
 

-The following program solves a  nxn  linear system  A.X = B  where  A = [ ai,j ]  is a tridiagonal matrix  i-e  ai,j = 0  if  | i - j | > 1
-In other words, only the main diagonal and its 2 nearest neighbors have non-zero elements.

-Gaussian elimination is used without pivoting, so the method can fail, for instance if  a11 = 0
  but practically, most of the problems leading to tridiagonal systems do not have these pitfalls.

-The components of the 3 diagonals are to be stored in column order into contiguous registers, starting with a register Rbb of your own choosing ( with bb  > 00 )
  followed by the  bi 's ( the n elements of B )

-The determinant  det A  is also computed and stored in R00.
 

Data Registers:    R00 = det A  at the end                                        ( Registers Rbb thru Ree are to be initialized before executing "3DLS" )
 

                             Rbb     = a11      Rbb+2 = a12                                                                                                                       R3n+bb-2 = Ree-n+1 = b1
                             Rbb+1 = a21     Rbb+3 = a22      Rbb+5 = a23                                                                                             R3n+bb-1 = Ree-n+2 = b2
                                                       Rbb+4 = a32      Rbb+6 = a33         ............
                                                                                  Rbb+7 = a43         ............                                                                                  .....................
                                                                                                                ............

                                                                                   R3n+bb-10 = an-3,n-2
                                                                                   R3n+bb-9 = an-2,n-2        R3n+bb-7 = an-2,n-1
                                                                                   R3n+bb-8 = an-1,n-2        R3n+bb-6 = an-1,n-1      R3n+bb-4 = an-1,n
                                                                                                                            R3n+bb-5 = an,n-1        R3n+bb-3 = an,n         R4n+bb-3 = Ree = bn
 

     >>>>  When the program stops, the solutions  x1 , .... , xn  have replaced  b1 , .... , bn  in registers  R3n+bb-2  thru  R4n+bb-3

Flags: /
Subroutines: /
 

   ( SIZE 4n+b-2 )
 
 
       STACK         INPUTS       OUTPUTS
            X     (bbb.eee)system    (bbb.eee)solution
            L              /              n

      (bbb.eee)system  is the control number of the system  ,  eee =  4n + bbb - 3  , bb > 00 ,  n = number of unknowns ,  n > 1
      (bbb.eee)solution is the control number of the solution

Example:                                                                                   With  bb = 01,  store                                            into

   2.x + 5.y                                     = 2                                       2   5                           2                 R01   R03                                           R17
   3.x + 7.y + 4.z                            = 4                                       3   7   4                      4                 R02   R04   R06                                  R18
               y + 3.z + 7.t                    = 7                                            1   3   7                 7                          R05   R07   R09                         R19
                     2.z + 4.t + 6.u           = 1                                                 2   4   6            1                                    R08   R10   R12                R20
                              8.t +  u   + 7.v  = 5                                                      8   1   7       5                                             R11   R13   R15       R21
                                      9.u + 4.v  = 6                                                           9   4       6                                                       R14   R16       R22

-Then,   1.022  XEQ "3DLS"  >>>>   17.022    ( in 9 seconds )

        R17 = x = -16.47810408            R20 = t = -0.480422463
        R18 = y =     6.991241632          R21 = u =  0.112313241
        R19 = z =     1.123905204          R22 = v =  1.247295209

-We have also  R00 = det A = 3882

Notes:

-The execution time is approximately proportional to n.
-With bb = 01, this program can solve a tridiagonal system of 80 equations in 80 unknowns in about  118 seconds.
-There is a risk of OUT OF RANGE  line 27 if the determinant exceeds 9.999999999 E99 - especially for large n-values. Set flag F24 in this case.

-"DIFF" calls "3DLS" as a subroutine.
 

One Dimensional Diffusion Equation ( Crank-Nicolson Method )
 

 "DIFF" solves the partial differential equations   T/x = A(x,y) 2T/y2 + B(x,y) 2T/xy + C(x,y) T/y + D(x,y) T + E(x,y)

     where  A , B , C , D , E  are known functions of  x and y

  with the initial values:               T(0,y) = F(y)
  and the boundary conditions:

               T(x,0)    = f(x)
             T(x,N.k)  = g(x)

-The interval  [0,N.k]  is divided into N parts.

    y

    |
    |----------------------------------------------------- N k
    |
    |
    |
 --|----------------------------------------------------- x
   0

-The partial derivatives are approximated by formulas of order 2 at the point (x+h/2,y)  ( Crank-Nicolson method )
 

                   T/x    ~  [ T(x+h,y) - T(x,y) ]/h
                 2T/xy ~ [ T(x+h,y+k) + T(x,y-k) - T(x+h,y-k) - T(x,y+k) ] / (2h.k)
 

-And the averages:

         T/y  ~ {  [ T(x,y+k) - T(x,y-k) ]/(2k)  +  [ T(x+h,y+k) - T(x+h,y-k) ]/(2k)  } / 2
        2T/y2  ~  {  [ 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

  the formulas are of order 2 both in space and time.

-  "DIFF"  must solve a tridiagonal linear system of M-1 equations at each step.
 

Data Registers:             R00 = x0                                      ( Registers R00 thru R05 are to be initialized before executing "DIFF" )

                                        R01 = h           R03 = k              R05 thru R13: temp     R14 .... = T(x,y)
                                        R02 = M         R04 = N > 2
 

       >>>>  When the program stops,   R14 = T(x,0)   R15 = T(x,k)   R16 = T(x,2k) .................  R14+M = T(x,N.k)          with  x = x0 + M.h

Flag:  F10 is cleared by this routine
Subroutines:  "3DLS" and:

   A program which computes A(x,y) assuming x is in X-register and t is in Y-register upon entry  named  LBL "AXY"
   A program which computes  B(x,y) assuming x is in X-register and t is in Y-register upon entry  named  LBL "BXY"
   A program which computes  C(x,y) assuming x is in X-register and t is in Y-register upon entry  named  LBL "CXY"
   A program which computes  D(x,y) assuming x is in X-register and t is in Y-register upon entry  named  LBL "DXY"
   A program which computes  E(x,y) assuming x is in X-register and t is in Y-register upon entry  named  LBL "EXY"

   A program which computes  T(x,0) assuming x is in X-register upon entry  named  LBL "TX0"
   A program which computes  T(x,L) assuming x is in X-register upon entry  named  LBL "TX1"
   A program which computes  T(0,y)  assuming y is in X-register upon entry  named  LBL "T0Y"

-We could have used data registers to store the names of these subroutines.
-Here, we use less registers.

   ( SIZE 5M+009 )
 
 
      STACK        INPUTS      OUTPUTS
           Y             /        bbb.eee
           X             /         x0+M.h

    bbb.eee = control number of the solution,      bbb = 14  ,   eee = 14+N

Example:     With the equation: T/x = y 2T/y2 - x 2T/xy + x y2 T/y + y2 T - 2 y2 exp ( - x y2 )  and

    x0 = 0                                    T(x,0) = 1
    h = 0.1   M = 10                    T(x,1) = exp(-x)
    k = 0.1   N = 10                    T(0,y) = 1

-Load the subroutines in main memory, for instance:
 
 
 01 LBL "AXY"
 02 X<>Y
 03 RTN
 04 LBL "BXY"
 05 CHS
 06 RTN
 07 LBL "CXY"
 08 X<>Y
 09 X^2
 10 *
 11 RTN
 12 LBL "DXY"
 13 X<>Y
 14 X^2
 15 RTN
 16 LBL "EXY"
 17 X<>Y
 18 X^2
 19 CHS
 20 STO Z       
 21 *
 22 E^X
 23 *
 24 ST+ X
 25 RTN
 26 LBL "TX0"
 27 CLX
 28 SIGN
 29 RTN
 30 LBL "TX1"
 31 CHS
 32 E^X
 33 RTN
 34 LBL "T0Y"
 35 CLX
 36 SIGN
 37 RTN
 38 END

 
    0   STO 00
  0.1  STO 01   1   STO 02
  0.1  STO 03  10  STO 04
 

   XEQ "DIFF"  >>>>    x = 0.1                                                 ---Execution time = 54s---
                         X<>Y   14.024

   and T(0.1,y)  are in registers  R14 to R24  ( y = 0 , 0.1 , ......... , 1 )

-Simply press R/S ( or XEQ 01 ) to continue with the same h and M ( or modify these parameters in R01 & R02 if needed ), it yields:
 
 
     x\y       0    1/10    2/10    3/10    4/10    5/10    6/10    7/10    8/10   9/10     L=1
      0       1       1       1       1       1       1       1       1       1       1       1
   1/10       1  0.9990  0.9960  0.9911  0.9842  0.9755  0.9649  0.9525  0.9384  0.9255  0.9048
   2/10       1  0.9980  0.9921  0.9823  0.9687  0.9515  0.9309  0.9070  0.8802  0.8506  0.8187
     ...    .......    .......    .......    .......    .......    .......    .......    .......   .........   .........    .......
      1       1  0.9900  0.9606  0.9136  0.8516  0.7781  0.6970  0.6120  0.5268  0.4446  0.3679
     reg.    R14    R15    R16    R17    R18    R19    R20    R21    R22    R23    R24

 
-The maximum error for x = 1  is about  -0.0007
 

Note:

-The exact solution is  T(x,y) = exp ( -x.y2 )
 

Two Dimensional Diffusion Equation
 

 "2DIFF" solves the partial differential equations

     T/x = A 2T/y2 + B 2T/yz + C 2T/z2 + D 2T/xy + E 2T/xz + F T/y + G T/z + H T + I

     where  A , B , C , D , E , F , G , H , I  are known functions of  x , y , z

  with the initial values:               T(0,y,z)
  and the boundary conditions:

             T(x,y,0)           T(x,0,z)
             T(x,y,L.l)         T(x,N.k,z)
 

Data Registers:             R00 = x0                                      ( Registers R00 thru R07 are to be initialized before executing "2DIFF" )

                                        R01 = h           R03 = k               R05 = l                 R07 thru R23: temp     R24 .... = T(x,y,z)
                                        R02 = M         R04 = N > 2        R06 = L > 2
 

Flag:  F00            CF 00 ->  relaxation parameter = LN(pi)
                             SF 00  -> relaxation parameter  =  1

Subroutines:     "AXYZ"  "BXYZ"  ..........................   "IXYZ"   that take x , y , z in registers X , Y , Z respectively and return the result in X-register
                           "T0YZ"   "TX0Z"  "TX1Z"  "TXY0"   "TXY1"   that take the 2 arguments ( x , y or x , z or y , z ) in registers X and Y

    ( SIZE 024+2 ( N+1) ( L+1 ) )
 
 
      STACK        INPUTS      OUTPUTS
           Y             /     24.eee.N+1
           X             /         x0+M.h

    24.eee.N+1 = control number of the solution

Example:      T/x =  2T/y2 + 2 2T/yz +  2T/z2 + y 2T/xy - z 2T/xz - y T/y + z T/z - ( y + z )2 T + exp( -x - y.z )

 with   T(0,y,z) = exp( -y.z )     T(x,0,z) = exp(-x)           T(x,y,0) = exp(-x)
                                               T(x,1,z) = exp(-x-z)         T(x,y,1) = exp(-x-y)
 

-Load these routines in main memory
 
 
 
 01 LBL "AXYZ"
 02 CLX
 03 SIGN
 04 RTN
 05 LBL "BXYZ"
 06 2
 07 RTN
 08 LBL "CXYZ"
 09 CLX
 10 SIGN
 11 RTN
 12 LBL "DXYZ"
 13 X<>Y
 14 RTN
 15 LBL "EXYZ"
 16 X<> Z
 17 CHS
 18 RTN
 19 LBL "FXYZ"
 20 X<>Y
 21 CHS
 22 RTN
 23 LBL "GXYZ"
 24 X<> Z
 25 RTN
 26 LBL "HXYZ"
 27 RDN
 28 +
 29 X^2
 30 CHS
 31 RTN
 32 LBL "IXYZ"
 33 X<> Z
 34 *
 35 +
 36 CHS
 37 E^X
 38 RTN
 39 LBL "T0YZ"
 40 *
 41 CHS
 42 E^X
 43 RTN
 44 LBL "TX0Z"
 45 CHS
 46 E^X
 47 RTN
 48 LBL "TX1Z"
 49 +
 50 CHS
 51 E^X
 52 RTN
 53 LBL "TXY0"
 54 CHS
 55 E^X
 56 RTN
 57 LBL "TXY1"
 58 +
 59 CHS
 60 E^X
 61 RTN
 62 END

 
    Then   0  STO 00  and if you choose  h = 0.1  M = 1 ,  k = 0.25  N = 4 ,  l = 0.2    L = 5

    0.1  STO 01         0.25  STO 03        0.2  STO 05
      1   STO 02           4     STO 04          5   STO 06

-The linear system is solved by a relaxation method, so the contents of many registers are taken as initial guesses
-If uou choose 0 as initial guesses,  SIZE 010   SIZE 084 or more
-The precision is controlled by the display format, let's choose FIX 4

  XEQ "2DIFF"  the sucessive variations are displayed and finally we get

                          x = 0.1   X<>Y   24.05305                                                   ---Execution = a few seconds with V41---

  And we have in registers  R24 to R53  T( 0.1 , y , z )

  y\z           0              1/5           2/5           3/5             4/5              1

  0        0.9048       0.9048     0.9048      0.9048      0.9048       0.9048          R24     R29     R34     R39     R44     R49
 1/4      0.9048       0.8603     0.8180      0.7778      0.7394       0.7047          R25     R30     R35     R40     R45     R50
 2/4      0.9048       0.8183     0.7401      0.6695      0.6053       0.5488    =    R26     R31     R36     R41     R46     R51
 3/4      0.9048       0.7784     0.6698      0.5764      0.4959       0.4274          R27    R32     R37     R42     R47     R52
   1       0.9048       0.7408     0.6065      0.4966      0.4066       0.3329          R28     R33     R38     R43     R48     R53

-The red values result from the differential equation
-The black values come from the boundary or initial conditions

-Press R/S or XEQ 10  to continue with the same h and M  ( or modify these values in R01 and R02 )
 

Notes:

-The exact solution is T(x,y,z) = exp( -x - y.z )
-The program is slow with a real HP41  ( about 16 minutes in FIX 3 )

-If you continue until x = 1 , you get T(1,y,z)  in the same registers ( R24 thru R53 )

  y\z           0              1/5           2/5           3/5             4/5              1

  0        0.3679       0.3679     0.3679      0.3679      0.3679       0.3679          R24     R29     R34     R39     R44     R49
 1/4      0.3679       0.3496     0.3324      0.3162      0.3010       0.2865          R25     R30     R35     R40     R45     R50
 2/4      0.3679       0.3325     0.3006      0.2718      0.2462       0.2231    =    R26     R31     R36     R41     R46     R51
 3/4      0.3679       0.3164     0.2721      0.2340      0.2014       0.1738          R27     R32     R37     R42     R47     R52
   1       0.3679       0.3012     0.2466      0.2019      0.1653       0.1353          R28     R33     R38     R43     R48     R53

-The maximum error is about  -0.0007
-Smaller stepsizes would give more accurate results.
 

Non-Linear Diffusion Equation ( Explicit Method )
 

"NLDIFF" try to solve the partial differential equation  T/x = f ( x , y , T , T/y , 2T/y2 )    where f  is a known function

  with the initial values:               T(0,y)

  and the boundary conditions:

             T(x,0)
             T(x,L)

-The interval  [0,N.k]  is divided into N parts.

    y

    |
    |----------------------------------------------------- N k
    |
    |
    |
 --|----------------------------------------------------- x
   0

-The partial derivatives are approximated by formulas of order 2, except T/x    ~  [ T(x+h,y) - T(x,y) ]/h   which is of order 1

-The method is often unstable and it's usually preferable to replace T(x,y) in the formula above by  [ T(x,y+k) + T(x,y-k) ] / 2   ( Lax method )

   This is applied if flag F00 is clear.
   But you can also use the 1st formula: simply set flag 00 ( SF 00 )
 

Data Registers:             R00 = x0                                     ( Registers R00 thru R05 are to be initialized before executing "NLDIFF" )

                                        R01 = h           R03 = k             R05 = y    R07 = T/y        R09 thru R11: temp     R12 .... = T(x,y)
                                        R02 = M         R04 = N            R06 = T   R08 = 2T/y2
 

       >>>>  When the program stops,   R12 = T(x,0)   R13 = T(x,k)   R14 = T(x,2k) .................  R12+M = T(x,N.k)          with  x = x0 + M.h

Flag:                     CF 00  Lax-Method
                              SF 00  No Lax-method ( more risky )

Subroutine:  A program,  LBL "dT/dX" that calculates  T/x = f ( x , y , T , T/y , 2T/y2 )
                                                               with  R00 = x , R05 to R08 = y , T , T/y , 2T/y2 respectively

     ( SIZE 014+2.N )
 
 
      STACK        INPUTS      OUTPUTS
           Y             /        12.eee
           X             /        x0+M.h

    12.eee = control number of the solution,      eee = 12+N

Example:    T/x = T ( - y + x T/y + 2T/y2 )       x0 = 0     T(0,y) = 1 = T(x,0)   T(x,1) = exp(-x)
 

-Load this routine in main memory:
 
 
 01 LBL "dT/dX"
 02 RCL 00
 03 RCL 07
 04 *
 05 RCL 05
 06 -
 07 RCL 08       
 08 +
 09 RCL 06
 10 *
 11 RTN
 12 LBL "TX0"  
 13 CLX
 14 SIGN
 15 RTN
 16 LBL "T0Y"  
 17 CLX
 18 SIGN
 19 RTN
 20 LBL "TX1"
 21 CHS
 22 E^X
 23 RTN            
 24 END

 
-Assuming you choose  h = 0.02   k = 0.2    M = 5   N = 5

     0     STO 00
   0.02  STO 01
     5     STO 02    STO 04
    0.2   STO 03

   SF 00  XEQ "NLDIFF"  >>>>   0.1    X<>Y   12.017   ( the HP41 displays the successive x-values )                ---Execution time = 56s---

  And we find in R12 thru R17     1  0.980   0.961   0.942   0.923   0.905

-Press R/S to continue with the same h and M-values ( or modify them in R01 & R02 ) and you gradually get:

y

 0        1    1            1            1              .............................      1
0.2      1    0.980     0.961     0.942       .............................      0.8184
0.4      1    0.961     0.923     0.887       .............................      0.6699
0.6      1    0.942     0.887     0.835       .............................      0.5483
0.8      1    0.923     0.852     0.786       .............................      0.4490
  1       1    0.905     0.819     0.741       .............................      0.3679

x =      0      0.1         0.2         0.3         .............................         1

-The exeact solution is  T(x,y) = exp( -x.y )  and the maximum error ( for x = 1 ) is about  -0.0005

-Surprisingly, Lax-method ( CF 00 ) is less accurate... and even unstable in this example.
 

Non-Linear 2D Diffusion Equation ( Explicit Method )
 

-The same method is used to solve   T/x = f ( x , y , z , T , T/y , T/z ,  2T/y2 , 2T/yz , 2T/z2  )

  knowing   T(0,y,z)   T(x,y,0)   T(x,0,z)  T(x,y,1)  T(x,1,z)        [  0 may be x0  and so on... ]
 
 

Data Registers:             R00 = x0                                      ( Registers R00 thru R07 are to be initialized before executing "NL2DIF" )

                                        R01 = h           R03 = k          R05 = l        R07 = y     R09 = T          R11 = T/z       R13 = 2T/yz
                                        R02 = M         R04 = N         R06 = L       R08 = z     R10 = T/y   R12 = 2T/y2    R14 = 2T/z2

                                         R15 thru R20: temp     R21 .... = T(x,y,z)
 

Flag:  F00            CF 00 ->  Lax-method
                             SF 00  -> No Lax-method

Subroutines:    "dT/dX"   that calculates  T/x  in X-register from  x , y , z , T , T/y , T/z ,  2T/y2 , 2T/yz , 2T/z2 in R00 R07 ... R14 respectively
                           "T0YZ"   "TX0Z"  "TX1Z"  "TXY0"   "TXY1"   that take the 2 arguments ( x , y or x , z or y , z ) in registers X and Y

    ( SIZE 021+2 ( N+1) ( L+1 ) )
 
 
      STACK        INPUTS      OUTPUTS
           Y             /     21.eee.N+1
           X             /        x0+M.h

    21.eee.N+1 = control number of the solution

Example:      T/x = T ( -y.z + 2 y T/y - 4 z T/z - y2 2T/y2 + 2 y.z 2T/yz - z2 2T/z2

  such that  x0 = 0    T(0,y,z) = 1

       T(x,y,0) = 1 = T(x,0,z)
       T(x,y,1) = exp( -x.y )
       T(x,1,z) = exp( -x.z )

-Load these subroutines in main memory:
 
 
 01 LBL "dT/dX"
 02 RCL 10
 03 RCL 07
 04 *
 05 ST+ X
 06 RCL 07
 07 RCL 08
 08 *
 09 -
 10 RCL 08
 11 RCL 11
 12 *
 13 4
 14 *
 15 -
 16 RCL 12
 17 RCL 07        
 18 X^2
 19 *
 20 -
 21 RCL 13
 22 RCL 07
 23 RCL 08
 24 *
 25 *
 26 ST+ X
 27 +
 28 RCL 14
 29 RCL 08
 30 X^2
 31 *
 32 -
 33 RCL 09
 34 *
 35 RTN
 36 LBL "T0YZ"
 37 CLX
 38 SIGN
 39 RTN
 40 LBL "TX0Z"
 41 CLX
 42 SIGN
 43 RTN
 44 LBL "TXY0"
 45 CLX
 46 SIGN
 47 RTN
 48 LBL "TX1Z"
 49 *
 50 CHS
 51 E^X
 52 RTN
 53 LBL "TXY1"
 54 *
 55 CHS
 56 E^X
 57 RTN
 58 END

 
-With  h = 0.01  M = 100  k = 0.2  N = 5  l = 0.2  L = 5

   0     STO 00
 0.01  STO 01  100  STO 02
  0.2   STO 03  STO 05
   5     STO 04  STO 06

  CF 00   XEQ "NL2DIF"  >>>>   1  X<>Y  21.05606                     ---Several seconds with V41 TURBO mode---

-You get T(1,y,z)  in registers  R21 thru R56

  y\z           0          1/5           2/5           3/5             4/5             1

  0             1           1              1               1               1               1              R21     R27     R33     R39     R45     R51
 1/5           1       0.969       0.939        0.909        0.877         0.819          R22     R28     R34     R40     R46     R52
 2/5           1       0.938       0.878        0.821        0.762         0.670    =    R23     R29     R35     R41     R47     R53
 3/5           1       0.906       0.817        0.736        0.657         0.549          R24     R30     R36     R42     R48     R54
 4/5           1       0.869       0.751        0.649        0.557         0.449          R25     R31     R37     R43     R49     R55
  1             1       0.819       0.670        0.549        0.449         0.368          R26     R32     R38     R44     R50     R56

-The exact solution is  T(x,y,z) = exp( -x.y.z )

Note:

-If you set F00 ( no Lax method ), you get OUT OF RANGE far before reaching x = 1  ( unstability )
-Even with Lax-method, the precision is nt very good:

    for example, T(1,2/5,3/5) ~ 0.821  whereas the correct result is about 0.787
    so the error is larger than 0.03
 

2D Poisson Equation
 

 "2DPOIS" solves     A 2U/x2 + B 2UT/xy + C 2U/y2 = D U/x + E U/y + F U + G

  with the boundary conditions   U(x,0)  U(0,y)   U(x,M.h)  U(N.k,y)
 

-We have to solve a linear system and this program uses a relaxation method.
-The precision of the solution ( of the linear system ) is controlled by the display format.
 
 

Data Registers:              R00 = temp                                  ( Registers R01 thru R04 are to be initialized before executing "2DPOIS" )

                                        R01 = h           R03 = k              R05 thru R15: temp     R16 .... = U(x,y)
                                        R02 = M         R04 = N
 

Flag:  F00             CF 00 ->  Relaxation parameter = Ln(PI)
                              SF 00 ->  Relaxation parameter =  1

Subroutines:   "AXY"  "BXY"  "CXY"  "DXY"  "EXY"  "FXY"  "GXY"   which take  x and y  in registers X and Y respectively
                          and  "UX0"  "U0Y"  "UX1"  "U1Y"

      ( SIZE 016+(M+1).(N+1) )
 
 
      STACK        INPUTS      OUTPUTS
           X             /     16.eee.N+1

    21.eee.N+1 = control number of the solution

Example:       2U/x2 + 2UT/xy + 2U/y2 = -y U/x - x U/y + x.y U - exp(-x.y)

  with  U(x,0) = 1 = U(0,y)     U(x,1) = exp(-x)    U(1,y) = exp(-y)
 

-Load these subroutines:
 
 
 01 LBL "AXY"
 02 CLX
 03 SIGN
 04 RTN
 05 LBL "BXY"
 06 CLX
 07 SIGN
 08 RTN
 09 LBL "CXY"
 10 CLX
 11 SIGN
 12 RTN
 13 LBL "DXY"
 14 X<>Y
 15 CHS
 16 RTN
 17 LBL "EXY"
 18 CHS
 19 RTN
 20 LBL "FXY"
 21 *
 22 RTN
 23 LBL "GXY"
 24 *
 25 CHS
 26 E^X
 27 CHS
 28 RTN
 29 LBL "UX0"
 30 CLX
 31 SIGN
 32 RTN
 33 LBL "U0Y"
 34 CLX
 35 SIGN
 36 RTN
 37 LBL "UX1"
 38 CHS
 39 E^X
 40 RTN
 41 LBL "U1Y"
 42 CHS
 43 E^X
 44 RTN
 45 END

 
-If you choose  h = 0.25  M = 4     k = 0.2   N = 5

  0.25  STO 01   4  STO 02
   0.2  STO 03    5  STO 04

  FIX 4  XEQ "2DPOIS"  >>>>  ( the HP41 displays the successive maximum corrections in the registers ) and finally  16.04505

  and you get in R16 thru R45

  x\y           0              1/5           2/5           3/5             4/5              1

  0             1                1              1               1               1                1              R16     R21     R26     R31     R36     R41
 1/4           1           0.9517     0.9052      0.8608      0.8186       0.7788          R17     R22     R27     R32     R37     R42
 2/4           1           0.9051     0.8189      0.7407      0.6701       0.6065    =    R18     R23     R28     R33     R38     R43
 3/4           1           0.8608     0.7408      0.6374      0.5485       0.4724           R19     R24     R29     R34     R39     R44
   1            1           0.8187     0.6703      0.5488      0.4493       0.3679           R20     R25     R30     R35     R40     R45

-The red values result from the differential equation
-The black values come from the boundary conditions

-The exact solution is  U(x,y) = exp(-x.y)
-So the errors remain smaller than 0.0005
 

Wave Equation ( Explicit Method )
 

-"WAVE"  solves equations of the type:     ¶2W/x2 = A(x,y) 2W/y2 + B(x,y) W/x + C(x,y) W/y + D(x,y) W + E(x,y)

  with the boundary conditions:       W(0,y)    W/x(0,y)
 

Data Registers:             R00 = x0                                     ( Registers R00 thru R04 are to be initialized before executing "WAVE" )

                                        R01 = h           R03 = k              R05 thru R13: temp     R14 .... = W(x,y)
                                        R02 = M         R04 = N             M <= INT(N/2)
 

Flags: /
Subroutines:   "AXY"  "BXY"  "CXY"  "DXY"  "EXY"  which take x & y in registers X & Y respectively.
                          "W0Y" & "dW0Y"

     ( SIZE 015+3.M )
 
 
      STACK        INPUTS      OUTPUTS
           T             /    bbb.eee(s-1)
           Z             /     bbb.eee(s)
           Y             /        ymin(s)
           X             /        x0+M.h

     bbb.eee(s) = control number of the solution at step s
  bbb.eee(s-1) = control number of the solution at step s-1

Example:     2W/x22W/y2 - y W/x + y W/y - x2 W + x.y exp(-x.y)

    with  W(0,y) = 1  and  W/x(0,y) = - y      (  x0 = 0 )

-Load the subroutines in main memory:
 
 
 01 LBL "AXY"
 02 CLX
 03 SIGN
 04 RTN
 05 LBL "BXY"
 06 X<>Y
 07 CHS
 08 RTN
 09 LBL "CXY"
 10 X<>Y
 11 RTN
 12 LBL "DXY"
 13 X^2
 14 CHS
 15 RTN
 16 LBL "EXY"
 17 *
 18 ENTER
 19 CHS
 20 E^X
 21 *
 22 RTN
 23 LBL "W0Y"
 24 CLX
 25 SIGN
 26 RTN
 27 LBL "dW0Y"
 28 CHS
 29 RTN
 30 END

 
-With h = 0.1 , M = 2  and  k = 0.1 ,  N = 10

    0    STO 00
  0.1   STO 01  STO 03
    2    STO 02
    10  STO 04  XEQ "WAVE"  >>>>        0.2
                                                 RDN         0.2
                                                 RDN      27.033
                                                 RDN      15.023

 And you get in R15 thru R23  and  in R27 to R33

y\x             0.1                      0.2                                   0.3                      0.4                                 0.5

 .9     R15 = 0.9900
 .8     R16 = 0.9800    R27 = 0.9603
 .7     R17 = 0.9700    R28 = 0.9408                    R17 = 0.9122
 .6     R18 = 0.9600    R29 = 0.9215                    R18 = 0.8842     R29 = 0.8481
 .5     R19 = 0.9500    R30 = 0.9023                    R19 = 0.8568     R30 = 0.8131                 R30 = 0.7712
 .4     R20 = 0.9400    R31 = 0.8834                    R20 = 0.8298     R31 = 0.7790
 .3     R21 = 0.9300    R32 = 0.8646                    R21 = 0.8034
 .2     R22 = 0.9200    R33 = 0.8460
 .1     R23 = 0.9100

-The column on the left ( for x = 0.1 ) is evaluated without the differential equation: only with W(0,y) & W/x(0,y)
-In the other columns, the differential equation plays the main role.
-We always have to save 2 columns to calculate the next step.

-Press R/S with the same parameters and it yields   0.4
                                                                   RDN  0.4
                                                                   RDN  29.031
                                                                   RDN  17.021

-The result are given in the next 2 columns above

-Then, we can only choose M = 1 to get the last result on the right!

    1  STO 02   R/S   >>>>   0.5
                                RDN    0.5
                                RDN    30.030
                                RDN    18.020   ( the previous results in R29-R30-R31 are now stored in R18-R19-R20 )

-So, we have estimated  W(0.5,0.5) ~ 0.7712

-The exact solution is  W(x,y) = exp(-x.y)   thus,  W(0.5,0.5) = 0.7788  and the error of our estimated value is about -0.0076

Notes:

-With x0 = 0   h = 0.05  M =10   k = 0.05  N = 20   you'll get  W(0.5,0.5) ~ 0.7750
-With x0 = 0   h = 0.02  M =25   k = 0.02  N = 50   you'll get  W(0.5,0.5) ~ 0.7773  ( in R90 )

  which suggests convergence to the exact result.

Warning:

-The formula will probably be unstable if B > 0
 

2D Wave Equation ( Explicit Method )
 

-The same method as above is used by "2DWAVE" to solve the equation:

   2W/x2 = A(x,y,z) 2W/y2 + B(x,y,z) 2W/yz + C(x,y,z) 2W/z2 + D(x,y,z) W/x + E(x,y,z) W/y + F(x,y,z) W/z + G(x,y,z) W + H(x,y,z)

   with the conditions    W(0,y,z)  and W/x(0,y,z)
 

Data Registers:            R00 = x0                                     ( Registers R00 thru R06 are to be initialized before executing "2DWAVE" )

                                        R01 = h           R03 = k           R05 = l          R07 thru R20: temp                          R21 .... = W(x,y,z)
                                        R02 = M         R04 = N          R06 = L        M <= Inf ( INT(N/2) , INT(L/2) )
 

Flags: /
Subroutines:   "AXYZ"  "BXYZ"  "CXYZ"  "DXYZ"  "EXYZ"   "FXYZ"  "GXYZ"  "HXYZ" which take x y z in registers X Y Z  respectively.
                          "W0YZ" & "dW0YZ"

     ( SIZE 032+3.(N.L-N-L) )
 
 
      STACK        INPUTS      OUTPUTS
           T             /    bbb.eee.ii(s)
           Z             /        zmin(s)
           Y             /        ymin(s)
           X             /        x0+M.h
           L             /   bbb.eee.ii(s-1)

     bbb.eee.ii(s) = control number of the solution at step s
  bbb.eee.ii(s-1) = control number of the solution at step s-1

Example:    2W/x2 = 2W/y2 + 2 2W/yz + 2W/z2 - W/x - y W/y + z W/z - ( y + z )2 W + 3 exp(-x-y.z)

    x0 = 0    W(0,y,z) = exp(-y.z) = - W/x(0,y,z)

-Load these subroutines:
 
 
 01 LBL "AXYZ"
 02 CLX
 03 SIGN
 04 RTN
 05 LBL "BXYZ"
 06 2
 07 RTN
 08 LBL "CXYZ"
 09 CLX
 10 SIGN
 11 RTN
 12 LBL "DXYZ"
 13 CLX
 14 SIGN
 15 CHS
 16 RTN
 17 LBL "EXYZ"
 18 X<>Y
 19 CHS
 20 RTN
 21 LBL "FXYZ"
 22 X<> Z
 23 RTN
 24 LBL "GXYZ"
 25 RDN
 26 +
 27 X^2
 28 CHS
 29 RTN
 30 LBL "HXYZ"
 31 X<> Z
 32 *
 33 +
 34 CHS
 35 E^X
 36 3
 37 *
 38 RTN
 39 LBL "W0YZ"
 40 *
 41 CHS
 42 E^X
 43 RTN
 44 LBL "dW0YZ"
 45 *
 46 CHS
 47 E^X
 48 CHS
 49 RTN
 50 END

 
-With h = 0.1  M = 2    k = 1/6  N = 6    l = 1/7  L = 7

    0   STO 00
  0.1  STO 01  2  STO 02
  1/6  STO 03  6  STO 04
   1/7 STO 05  7  STO 06

   XEQ "2DWAVE"   >>>>     0.2
                                 RDN      0.3333        ( in fact 2/6 )
                                 RDN      0.2857        ( in fact 2/7 )
                                 RDN      51.06203
                               LASTX    21.05005

  and we find  W( 0.2 , y , z )  in registers  R51 to R62
 

    y\z      2/7       3/7       4/7       5/7

   2/6    0.738   0.703   0.669   0.637         R51    R54    R57    R60
   3/6    0.706   0.657   0.611   0.568   =    R52    R55    R58    R61
   4/6    0.676   0.615   0.559   0.508         R53    R56    R59    R62

-If we continue with M = 1  STO 02

      R/S   >>>>   0.3   RDN   0.5   RDN   0.4286    RDN    33.03401    LASTX   21.03203

   And we get   W( 0.3 , y , z )  in registers  R33 & R34

    y\z      3/7      4/7

    0.5   0.601   0.558            ( The exact values are  0.5979  and  0.5567  respectively rounded to 4D )

-The exact solution is Wx,y,z) = exp( -x-y.z )

Note:

-If you store M = 3 in R02 at the beginning, execution time = 3m32s

Warning:

-The formula will probably be unstable if D > 0
 

Garden-Hose Method
 

 "SHOOT" solves the following boundary problem:

  Given the ordinary differential equation  y'' = T ( x , y , y' )
  Find the solution satisfying  y(a) = A & y(b) = B

  You give 2 initial guesses for y'(a)  which lead to 2 results for y(b)
  and your HP41 uses the secant method to find the correct y'(a) that will give y(b) = B

  This type of problem often arises when we search the minimum - or extremum - of the integral   I = §ab  L( x , y , y' ) dx  with y(a) = A & y(b) = B
  Your HP41 will also calculate this integral.

-You have, however, to do some calculus to get T ( x , y , y' )  from  L( x , y , y' )    ( Euler-Lagrange equations )

  "SHOOT"  uses the "classical" Runge-Kutta method of order 4 ( RK4 ) and Simpson's rule to evaluate the integral.
 

>>>  You have to store an even integer N in R00 ,

         h = (b-a) / N  will be the stepsize to integrate y'' = T ( x , y , y' ) and to evaluate I = §ab  L( x , y , y' ) dx

>>>  The precision is controlled by the display format.
 
 

Data Registers:            R00 = N  ( even integer )                        ( Registers R00 thru R04 are to be initialized before executing "SHOOT" )

                                        R01 = a           R03 = y(a)         R05 = y' (a)          R07 = I           R09 thru R19: temp         R20 .... = y(a) , y(a+h) , ..............
                                        R02 = b           R04 = y(b)         R06 = y' (b)         R08 = 20.eee
 

Flag:  F00    CF 00  the (N+1) y-values are stored in R20 R21 ........ R20+N
                     SF 00   the computed y-values are not stored

Subroutines:   LBL "T"  &  LBL "L"  to compute T ( x , y , y' )  & L ( x , y , y' )   from  x , y , y'  in registers X , Y , Z respectively.
 
 
 
      STACK        INPUTS      OUTPUTS
           T             /        20.eee
           Z             /          Imin
           Y             /         y' (b)
           X             /         y' (a)
           L             /       y(b) - B

     y(b) - B   should be a small number

Example:    Find the minimum of the integral  I = §-1+1  y ( 1 + y' 2 ) 1/2  dx    with  y(-1) = y(+1) = cosh 1 = 1.543080635....

-After some calculus, the Euler-Lagrange equation  L/y = ( d/dx ) ( L/y' )  gives

    y'' = ( 1 + y' 2 ) / y  =  T ( x , y , y' )
 
 
 01 LBL "T"
 02 X<> Z
 03 X^2
 04 1
 05 +
 06 X<>Y
 07 /
 08 RTN
 09 LBL "L"
 10 X<> Z
 11 X^2
 12 1
 13 +
 14  SQRT
 15 *
 16 RTN
 17 END

 
-If you choose N = 20    20  STO 00

   -1  STO 01   1   STO 02
    1.543080635    STO 03   STO 04

-If you choose  FIX 4  and  -0.5 & -0.6  as initial guesses for y'(a)

   FIX 4   CF 00
    -0.5  ENTER^
    -0.6  XEQ "SHOOT"  >>>>  -1.1752         ( exact result = - Sinh 1 )
                                       RDN     1.1752
                                       RDN     2.813445
                                       RDN     20.040
                                    LASTX    -2 E-9

-So, the minimum for the integral  I = §-1+1  y ( 1 + y' 2 ) 1/2  dx  [  if  y(-1) = y(+1) ]    is  2.813445

  y'(-1) = -1.1752  &  y'(1) = +1.1752

-And we find in R20 to R40 the following values for y(x):

   x   |     -1        -0.9       -0.8       -0.7      -0.6      -0.5      -0.4       -0.3      -0.2        -0.1         0        +0.1     +0.2      +0.3      +0.4     +0.5
--------------------------------------------------------------------------------------------------------------------------------------------------
   y   |  1.5431  1.4331  1.3374  1.2552  1.1855  1.1276  1.0811  1.0453  1.0201   1.0050  1.0000  1.0050  1.0201  1.0453  1.0811  1.1276
 

   x   |      +0.6      +0.7       +0.8      +0.9        +1
------------------------------------------------------
   y  |     1.1855   1.2552   1.3374   1.4331   1.5431
 

-The exact solution is  y(x) = Cosh x
-In fact, almost 7 decimals are correct.

-And the integral = 2.813430204  ( rounded to 9 decimals )
  so our approximation is correct to almost 5 decimals

-This curve minimizes the area of a surface of revolution ( multiply I by 2 PI )

Notes:

-If you store an odd number in R00, your HP41 will add 1 to this register.

-If you only want to solve y'' = T ( x , y , y' ) , delete lines 09 to 16 in the subroutine and add LBL "L" after line 01.
-Thus you will have another approximation of   §ab  T( x , y , y' ) = y'(b) - y'(a)
 

Euler-Lagrange Equation
 

"EULAG" finds the minimum - or extremum - of the integral   I = §ab  L( x , y , y' ) dx  with y(a) = A & y(b) = B
 without asking you for any calculus

-The Euler-Lagrange equation may be  written   L/y = 2L/xy' + y'  2L/yy' + y'' 2L/y'2

-So,  y'' may be expressed as a function of  x , y , y'  provided  2L/y'2 # 0

-The derivatives are approximated by formulas of order 2 and then, the HP41 evaluates

     y(x+h) ~ y(x) + h.y'(x) + ( h2 / 2 ) y''(x)
     y'(x+h) ~ y'(x) + h.y''(x)

-Like the program above, "EULAG" uses Simpson's rule to approximate the integral.
 

>>>  The precision is controlled by the display format.
 

Data Registers:            R00 = N  ( even integer )                        ( Registers R00 thru R04 are to be initialized before executing "EULAG" )

                                        R01 = a           R03 = y(a)         R05 = y' (a)          R07 = I           R09 thru R19: temp         R20 .... = y(a) , y(a+h) , ..............
                                        R02 = b           R04 = y(b)         R06 = y' (b)         R08 = 20.eee
 

Flag:  F00    CF 00  the (N+1) y-values are stored in R20 R21 ........ R20+N
                     SF 00   the computed y-values are not stored

Subroutine:    LBL "L"  to compute  L ( x , y , y' )   from  x , y , y'  in registers X , Y , Z  respectively.
 
 
 
      STACK        INPUTS      OUTPUTS
           T             /        20.eee
           Z             /          Imin
           Y             /         y' (b)
           X             /         y' (a)
           L             /       y(b) - B

     y(b) - B   should be a small number

Example:    The same one:  Find the minimum of the integral  I = §-1+1  y ( 1 + y' 2 ) 1/2  dx    with  y(-1) = y(+1) = cosh 1 = 1.543080635....

-Load this subroutine in main memory:
 
 
 01 LBL "L"
 02 X<> Z
 03 X^2
 04 1
 05 +
 06  SQRT
 07 *
 08 RTN
 09 END

 
-If you choose N = 20    20  STO 00

   -1  STO 01   1   STO 02
    1.543080635    STO 03   STO 04

-If you choose  FIX 4  and  -0.5 & -0.6  as initial guesses for y'(a)

   FIX 4   CF 00
    -0.5  ENTER^
    -0.6  XEQ "EULAG"  >>>>   -1.2157           ( exact result = -Sinh 1 = -1.175201194.... )
                                       RDN     1.1659           ( exact result = Sinh 1 = 1.175201194.... )
                                       RDN     2.8128
                                       RDN     20.040
                                    LASTX    1.3  E-07

-So, the minimum for the integral  I = §-1+1  y ( 1 + y' 2 ) 1/2  dx  [  if  y(-1) = y(+1) ]    is about  2.8128

  y'(-1) = -1.2157  &  y'(1) = +1.1659

-And we find in R20 to R40 the following values for y(x):

   x   |     -1        -0.9     -0.8     -0.7    -0.6    -0.5     -0.4     -0.3    -0.2     -0.1       0      +0.1    +0.2    +0.3    +0.4   +0.5
------------------------------------------------------------------------------------------------------------------------------
   y   |  1.5431  1.430   1.331  1.247  1.177  1.119  1.072  1.037  1.012   0.998  0.993  0.999  1.016  1.042  1.079  1.126
 

   x   |     +0.6     +0.7    +0.8     +0.9      +1
------------------------------------------------
   y  |     1.185   1.255   1.338   1.434   1.5431
 

-The results are of course less accurate
-The exact integral = 2.813430204  ( rounded to 9 decimals )  so, the error is about -0.0006  ( we've found 2.8128 )

-However, with  N = 100  STO 00  ( and V41 in turbo mode )  we get  I = 2.81341

Notes:

-If you store an odd number in R00, your HP41 will add 1 to this register.

-Set Flag 00 ( SF 00 ) if you don't want to store the y-values - or if there are not enough registers...
 
 

References:

[1]  Press , Vetterling , Teukolsky , Flannery - "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4
[2]  Evans , Blackledge , Yardley - Numerical Methods for Partial Differential Equations" - Springer - ISBN 3-540-76125-X
[3]  F. Scheid - "Numerical Analysis" - Mc Graw Hill   ISBN 07-055197-9