Partial Differential Equations for the HP41

Overview
 

 1°)  Tridiagonal Linear Systems
 2°)  One Dimensional Diffusion Equation ( Crank-Nicolson Method )
 3°)  Two Dimensional Diffusion Equation ( Crank-Nicolson Method )
 4°)  Non-Linear Diffusion Equation ( Explicit Method )
 5°)  2-Dimensional Non-Linear Diffusion Equation ( Explicit Method )
 6°)  2D Poisson Equation
 7°)  Wave Equation ( Explicit Method )
 8°)  2D Wave Equation ( Explicit Method )
 9°)  Garden-Hose Method
10°) Euler-Lagrange Equation
 

-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
   ¶²f/¶x² ~ [ f(x+h) - 2 f(x) + f(x+h) ] / h²

  ¶²f/¶x¶y ~ [ 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

>>>  All these programs are in a PARDIFFEQ module which may be downloaded  here
 

1°)  Tridiagonal Linear Systems
 

-The following program solves a  nxn  linear system  A.X = B  where  A = [ a_i,j ]  is a tridiagonal matrix  i-e  a_i,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  a₁₁ = 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  b_i '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     = a₁₁    •  Rbb+2 = a₁₂                                                                                                                     •  R3n+bb-2 = Ree-n+1 = b₁
                           •  Rbb+1 = a₂₁   •  Rbb+3 = a₂₂    •  Rbb+5 = a₂₃                                                                                           •  R3n+bb-1 = Ree-n+2 = b₂
                                                     •  Rbb+4 = a₃₂    •  Rbb+6 = a₃₃         ............
                                                                                •  Rbb+7 = a₄₃         ............                                                                                  .....................
                                                                                                                ............

                                                                                 •  R3n+bb-10 = a_n-3,n-2
                                                                                 •  R3n+bb-9 = a_n-2,n-2      •  R3n+bb-7 = a_n-2,n-1
                                                                                 •  R3n+bb-8 = a_n-1,n-2      •  R3n+bb-6 = a_n-1,n-1    •  R3n+bb-4 = a_n-1,n
                                                                                                                          •  R3n+bb-5 = a_n,n-1      •  R3n+bb-3 = a_n,n       •  R4n+bb-3 = Ree = b_n
 

     >>>>  When the program stops, the solutions  x₁ , .... , x_n  have replaced  b₁ , .... , b_n  in registers  R3n+bb-2  thru  R4n+bb-3

Flags: /
Subroutines: /
 
 

01  LBL "3DLS" 02  INT 03  LASTX 04  FRC 05   E3 06  * 07  RCL X 08  3 09  + 10  R^ 11  - 12  4 13  / 14  STO M 15  - 16  1

17  STO 00 18  + 19  STO N 20   E3 21  / 22  + 23  STO O 24  LBL 01 25  RCL O 26  RCL IND X 27  ST* 00 28  ST/ IND N 29  ISG O 30  ISG O 31  FS? 30 32  GTO 02

33  ST/ IND O 34  CLX 35  SIGN 36  + 37  RCL IND N 38  RCL IND Y 39  * 40  ISG N 41  CLX 42  ST- IND N 43  RCL IND O 44  LASTX 45  * 46  ISG O 47  ST- IND O 48  GTO 01

49  LBL 02 50  CLX 51  SIGN 52  - 53  INT 54  3 E-5 55  + 56  RCL M 57  DSE X 58  STO O 59  X<>Y 60  LBL 03 61  RCL IND X 62  RCL IND N 63  * 64  DSE N

65  ST- IND N  66  DSE Y 67  X<>Y 68  DSE O 69  GTO 03 70  RCL N 71  RCL M 72  + 73  DSE X 74   E3 75  ST/ Y 76  RDN 77  ST+ N 78  X<> N 79  CLA 80  END

 
   ( 132bytes / 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.
 

2°) One Dimensional Diffusion Equation ( Crank-Nicolson Method )
 

 "DIFF" solves the partial differential equations   ¶T/¶x = A(x,y) ¶²T/¶y² + B(x,y) ¶²T/¶x¶y + 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
               ¶²T/¶x¶y ~ [ 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
        ¶²T/¶y²  ~  {  [ T(x+h,t) - 2.T(x,t) + T(x+h,t) ]/h²  +  [ T(x+h,t+k) - 2.T(x,t+k) + T(x+h,t+k) ]/h² } /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 = x₀                                     ( 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 = x₀ + 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.
 
 

01 LBL "DIFF"  02 RCL 04  03 STO 09  04 LBL 00  05 RCL 03  06 RCL 09  07 *  08 XEQ "T0Y"       09 RCL 09  10 14  11 +  12 X<>Y  13 STO IND Y  14 DSE 09  15 GTO 00  16 CLX  17 XEQ "T0Y"  18 STO 14  19 LBL 01  20 RCL 02  21 STO 12  22 LBL 02  23 CLX  24 STO 08  25 RCL 04  26 4  27 *

28 10  29 +  30 STO 09  31  E3  32 /  33 15  34 STO 11  35 +  36 RCL 04  37 +  38 2 E-5  39 +  40 STO 10  41 RCL 01  42 ST+ 00  43 RCL 00  44 XEQ "TX0"       45 STO 05  46 SF 10  47 RCL 01  48 2  49 /  50 ST- 00  51 DSE 11  52 LBL 03  53 RCL 03  54 ST+ 08

55 RCL 08  56 RCL 00  57 XEQ "EXY"  58 ST+ X  59 CHS  60 STO IND 09  61 RCL 08  62 RCL 00  63 XEQ "CXY"  64 RCL 03  65 ST+ X  66 /  67 STO 07  68 RCL 08  69 RCL 00  70 XEQ "BXY"     71 RCL 01  72 RCL 03  73 *  74 /  75 STO 13  76 RCL 08  77 RCL 00  78 XEQ "AXY"  79 RCL 03  80 X^2  81 /

82 STO 06  83 RCL 07  84 -  85 STO Y  86 RCL 13  87 ST+ Z  88 -  89 FC?C 10  90 GTO 04  91 RCL 05  92 *  93 ST- IND 09  94 GTO 05  95 LBL 04  96 STO IND 10    97 ISG 10  98 LBL 05  99 X<>Y 100 RCL IND 11 101 * 102 ST- IND 09 103 RCL 08 104 RCL 00 105 XEQ "DXY" 106 RCL 06 107 ST+ 07 108 ST+ X

109 - 110 STO Y 111 RCL 01 112 1/X 113 ST+ X 114 ST+ Z 115 - 116 STO IND 10   117 X<>Y 118 ISG 11 119 CLX 120 RCL IND 11 121 * 122 ST- IND 09 123 ISG 11 124 CLX 125 RCL IND 11 126 RCL 07 127 RCL 13 128 ST+ 07 129 - 130 * 131 ST- IND 09 132 ISG 10 133 FS? 30 134 GTO 06 135 RCL 07

136 STO IND 10   137 CLX 138 SIGN 139 ST+ 09 140 ST- 10 141 ST- 11 142 GTO 03 143 LBL 06 144 RCL 01 145 2 146 / 147 ST+ 00 148 RCL 00 149 XEQ "TX1" 150 STO IND 11 151 RCL 07 152 * 153 ST- IND 09 154 RCL 00 155 X<> 05 156 STO 14 157 RCL 09 158  E3 159 / 160 15 161 + 162 RCL 04

163 + 164 XROM "3DLS" 165 RCL 05 166 STO 00 167 X<>Y 168 INT 169 .015 170 + 171 RCL 04 172 DSE X 173  E6 174 / 175 + 176 REGMOVE 177 VIEW 00 178 DSE 12 179 GTO 02 180 RCL 11 181  E3 182 / 183 14 184 + 185 RCL 00 186 RTN 187 GTO 01 188 END

 
      ( 304 bytes / 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 ¶²T/¶y² - x ¶²T/¶x¶y + x y² ¶T/¶y + y² T - 2 y² exp ( - x y² )  and

    x₀ = 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.y² )
 

3°)  Two Dimensional Diffusion Equation  ( Crank-Nicolson Method )
 

 "2DIFF" solves the partial differential equations

     ¶T/¶x = A ¶²T/¶y² + B ¶²T/¶y¶z + C ¶²T/¶z² + D ¶²T/¶x¶y + E ¶²T/¶x¶z + 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)

-The linear system is solved by a relaxation method.
-The initial content of these registers may change the rate of convergence.

-If you want to store 0 in the registers,  XEQ SIZE 007  before the correct SIZE...
 

Data Registers:           •  R00 = x₀                                     ( Registers R00 thru R06 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
 
 

01 LBL "2DIFF"   02 RCL 04  03 1  04 +  05 STO 09  06 RCL 06  07 LASTX  08 +  09 STO 10  10 *  11 STO 16  12 23  13 +  14 STO 13  15 RCL 03  16 RCL 04  17 *  18 STO 07  19 RCL 05  20 RCL 06  21 *  22 STO 08  23 LBL 00  24 RCL 08  25 RCL 07  26 XEQ "T0YZ"  27 STO IND 13  28 RCL 03  29 ST- 07  30 DSE 13  31 DSE 09  32 GTO 00  33 RCL 03  34 RCL 04  35 *  36 STO 07  37 RCL 05  38 ST- 08  39 RCL 04  40 1  41 +  42 STO 09  43 DSE 10  44 GTO 00  45 LBL 10  46 RCL 02  47 STO 17  48 LBL 01  49 RCL 04  50 1  51 +  52 RCL 06  53 ST* Y  54 +  55 .1  56 %  57 +  58  E3  59 /  60 24.024  61 +  62 REGMOVE  63 RCL 04  64 24.023  65 +  66 STO 15  67 RCL 03

68 RCL 04  69 *  70 STO 07  71 RCL 01  72 ST+ 00  73 LBL 04  74 RCL 07  75 RCL 00  76 XEQ "TXY0"  77 STO IND 15  78 RCL 03  79 ST- 07  80 DSE 15  81 GTO 04  82 RCL 16  83 23.023  84 +  85 RCL 04  86 1  87 +  88  E5  89 /  90 +  91 STO 14  92 RCL 04  93 -  94 STO 13  95 STO 15  96 RCL 05  97 RCL 06  98 *  99 STO 08 100 LBL 05 101 RCL 08 102 RCL 00 103 XEQ "TX0Z" 104 STO IND 15 105 RCL 05 106 ST- 08 107 DSE 15 108 GTO 05 109 RCL 03 110 RCL 04 111 * 112 STO 07 113 RCL 05 114 RCL 06 115 * 116 STO 08 117 LBL 06 118 RCL 08 119 RCL 00 120 XEQ "TX1Z" 121 STO IND 14 122 RCL 05 123 ST- 08 124 DSE 14 125 GTO 06 126 RCL 16 127 RCL 13 128 INT 129 DSE X 130  E3 131 / 132 + 133 23 134 +

135 STO 15 136 LBL 07 137 RCL 07 138 RCL 00 139 XEQ "TXY1" 140 STO IND 15 141 RCL 03 142 ST- 07 143 DSE 15 144 GTO 07 145 RCL 01 146 2 147 / 148 ST- 00 149 LBL 12 150 CLX 151 STO 22 152 RCL 04 153 RCL 06 154 ST* Y 155 + 156 DSE X 157 23 158 + 159 STO 21 160 RCL 06 161 DSE X 162 STO 10 163 RCL 05 164 * 165 STO 08 166 LBL 02 167 RCL 04 168 DSE X 169 STO 09 170 RCL 03 171 * 172 STO 07 173 LBL 03 174 RCL 08 175 RCL 07 176 RCL 00 177 XEQ "AXYZ" 178 RCL 03 179 X^2 180 ST+ X 181 / 182 STO 11 183 RCL 08 184 RCL 07 185 RCL 00 186 XEQ "BXYZ" 187 RCL 03 188 RCL 05 189 * 190 8 191 * 192 / 193 STO 12 194 RCL 08 195 RCL 07 196 RCL 00 197 XEQ "CXYZ" 198 RCL 05 199 X^2 200 ST+ X 201 /

202 STO 13 203 RCL 08 204 RCL 07 205 RCL 00 206 XEQ "DXYZ" 207 RCL 01 208 RCL 03 209 * 210 ST+ X 211 / 212 STO 14 213 RCL 08 214 RCL 07 215 RCL 00 216 XEQ "EXYZ" 217 RCL 01 218 RCL 05 219 * 220 ST+ X 221 / 222 STO 15 223 RCL 08 224 RCL 07 225 RCL 00 226 XEQ "FXYZ" 227 RCL 03 228 4 229 * 230 / 231 STO 18 232 RCL 08 233 RCL 07 234 RCL 00 235 XEQ "GXYZ" 236 RCL 05 237 4 238 * 239 / 240 STO 19 241 RCL 08 242 RCL 07 243 RCL 00 244 XEQ "HXYZ" 245 2 246 / 247 STO 20 248 RCL 21 249 1 250 + 251 RCL IND X 252 RCL 11 253 RCL 14 254 + 255 RCL 18 256 + 257 * 258 STO 23 259 RCL 21 260 DSE X 261 RCL IND X 262 RCL 11 263 RCL 14 264 - 265 RCL 18 266 - 267 * 268 ST+ 23

269 RCL 21 270 RCL 04 271 + 272 1 273 + 274 RCL IND X  275 RCL 13 276 RCL 15 277 + 278 RCL 19 279 + 280 * 281 ST+ 23 282 RCL 21 283 RCL 04 284 - 285 DSE X 286 RCL IND X 287 RCL 13 288 RCL 15 289 - 290 RCL 19 291 - 292 * 293 ST+ 23 294 RCL 21 295 RCL 04 296 2 297 + 298 + 299 RCL 21 300 LASTX 301 - 302 RCL 21 303 RCL 04 304 + 305 RCL 21 306 ST- L 307 RCL IND T 308 RCL IND T 309 RCL IND T 310 RCL IND L 311 + 312 - 313 + 314 RCL 12 315 * 316 ST+ 23 317 RCL 16 318 ST+ 21 319 RCL IND 21 320 RCL 01 321 1/X 322 RCL 11 323 RCL 13 324 + 325 ST+ X 326 - 327 RCL 20 328 + 329 * 330 ST+ 23 331 RCL 08 332 RCL 07 333 RCL 00 334 XEQ "IXYZ" 335 ST+ 23

336 RCL 21 337 1 338 + 339 RCL IND X 340 RCL 11 341 RCL 14 342 - 343 RCL 18 344 + 345 * 346 ST+ 23 347 RCL 21 348 DSE X 349 RCL IND X  350 RCL 11 351 RCL 14 352 + 353 RCL 18 354 - 355 * 356 ST+ 23 357 RCL 21 358 RCL 04 359 + 360 1 361 + 362 RCL IND X 363 RCL 13 364 RCL 15 365 - 366 RCL 19 367 + 368 * 369 ST+ 23 370 RCL 21 371 RCL 04 372 - 373 DSE X 374 RCL IND X 375 RCL 13 376 RCL 15 377 + 378 RCL 19 379 - 380 * 381 ST+ 23 382 RCL 21 383 RCL 04 384 2 385 + 386 + 387 RCL 21 388 LASTX 389 - 390 RCL 21 391 RCL 04 392 + 393 RCL 21 394 ST- L 395 RCL IND T 396 RCL IND T 397 RCL IND T 398 RCL IND L 399 + 400 - 401 + 402 RCL 12

403 * 404 ST+ 23 405 RCL 16 406 ST- 21 407 RCL 23 408 RCL 01 409 1/X 410 RCL 11 411 RCL 13 412 + 413 ST+ X 414 + 415 RCL 20 416 - 417 / 418 RCL IND 21 419 - 420 PI 421 LN 422 FS? 00 423 SIGN 424 * 425 ST+ IND 21 426 ABS 427 RCL 22 428 X<Y? 429 X<>Y 430 STO 22 431 RCL 03 432 ST- 07 433 DSE 21 434 DSE 09 435 GTO 03 436 RCL 05 437 ST- 08 438 DSE 21 439 DSE 21 440 DSE 10 441 GTO 02 442 RCL 22 443 RND 444 VIEW X 445 X#0? 446 GTO 12 447 RCL 01 448 2 449 / 450 ST+ 00 451 VIEW 00 452 DSE 17 453 GTO 01 454 24.023 455 RCL 16 456  E3 457 / 458 + 459 RCL 04 460 1 461 + 462  E5 463 / 464 + 465 RCL 00 466 RTN 467 GTO 10 468 END

 
    ( 719 bytes / 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 =  ¶²T/¶y² + 2 ¶²T/¶y¶z +  ¶²T/¶z² + y ¶²T/¶x¶y - z ¶²T/¶x¶z - y ¶T/¶y + z ¶T/¶z - ( y + z )² 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 007   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.
 

4°)  Non-Linear Diffusion Equation ( Explicit Method )
 

"NLDIFF" try to solve the partial differential equation  ¶T/¶x = f ( x , y , T , ¶T/¶y , ¶²T/¶y² )    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 = x₀                                    ( 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 = ¶²T/¶y²
 

       >>>>  When the program stops,   R12 = T(x,0)   R13 = T(x,k)   R14 = T(x,2k) .................  R12+M = T(x,N.k)          with  x = x₀ + 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 , ¶²T/¶y² )
                                                               with  R00 = x , R05 to R08 = y , T , ¶T/¶y , ¶²T/¶y² respectively
 
 

01 LBL "NLDIFF"  02 RCL 04  03 STO 09  04 LBL 00  05 RCL 03  06 RCL 09  07 *  08 XEQ "T0Y"  09 RCL 09  10 12  11 +  12 X<>Y  13 STO IND Y  14 DSE 09  15 GTO 00  16 CLX  17 XEQ "T0Y"  18 STO 12  19 LBL 01  20 RCL 02  21 STO 10  22 LBL 02  23 RCL 04

24 DSE X  25  E3  26 /  27 13.012  28 +  29 STO 09  30 CLX  31 STO 05  32 LBL 03  33 RCL 03  34 ST+ 05  35 RCL 09  36 RCL 09  37 1  38 ST+ Z  39 -  40 RCL IND Y      41 STO 08  42 RCL IND Y  43 ST+ 08  44 -  45 RCL 03  46 ST+ X

47 /  48 STO 07  49 RCL IND 09  50 STO 06  51 ST+ X  52 ST- 08  53 RCL 03  54 X^2  55 ST/ 08  56 XEQ "dT/dX"    57 RCL 01  58 *  59 RCL 09  60 RCL 09  61 1  62 ST+ Z  63 -  64 R^  65 RCL IND Z  66 RCL IND Z  67 +  68 2  69 /

70 FS? 00  71 RDN  72 FS? 00  73 RCL 06  74 +  75 RCL 09  76 RCL 04  77 +  78 1  79 +  80 X<>Y  81 STO IND Y      82 ISG 09  83 GTO 03  84 RCL 04  85 1  86 +  87  E6  88 /  89 1.012  90 +  91 RCL 09  92 INT

93 +  94 REGMOVE  95 RCL 01  96 ST+ 00  97 RCL 00  98 XEQ "TX0"  99 STO 12 100 RCL 00 101 XEQ "TX1" 102 STO IND 09   103 VIEW 00 104 DSE 10 105 GTO 02 106 RCL 09 107 INT 108  E3 109 / 110 12 111 + 112 RCL 00 113 RTN 114 GTO 01 115 END

 
       ( 193 bytes / 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 + ¶²T/¶y² )       x₀ = 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.
 

5°) 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 ,  ¶²T/¶y² , ¶²T/¶y¶z , ¶²T/¶z²  )

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

Data Registers:           •  R00 = x₀                                     ( 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 = ¶²T/¶y¶z
                                      •  R02 = M       •  R04 = N       •  R06 = L       R08 = z     R10 = ¶T/¶y   R12 = ¶²T/¶y²    R14 = ¶²T/¶z²

                                         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 ,  ¶²T/¶y² , ¶²T/¶y¶z , ¶²T/¶z² 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
 
 

01 LBL "NL2DIF"  02 RCL 04  03 1  04 +  05 STO 09  06 RCL 06  07 LASTX  08 +  09 STO 10  10 *  11 STO 20  12 20  13 +  14 STO 13  15 RCL 04  16 RCL 06  17 DSE X  18 ST* Y  19 +  20 RCL 04  21 +  22 STO 16  23 RCL 03  24 RCL 04  25 *  26 STO 07  27 RCL 05  28 RCL 06  29 *  30 STO 08  31 LBL 00  32 RCL 08  33 RCL 07  34 XEQ "T0YZ"  35 STO IND 13  36 RCL 03  37 ST- 07  38 DSE 13  39 DSE 09  40 GTO 00  41 RCL 03  42 RCL 04  43 *  44 STO 07  45 RCL 05

46 ST- 08  47 RCL 04  48 1  49 +  50 STO 09  51 DSE 10  52 GTO 00  53 LBL 10  54 RCL 02  55 STO 17  56 LBL 01  57 RCL 04  58 RCL 06             59 ST* Y  60 +  61 DSE X  62 20  63 +  64 STO 15  65 RCL 06  66 DSE X  67 STO 19  68 RCL 05  69 *  70 STO 08  71 LBL 02  72 RCL 04  73 DSE X  74 STO 18  75 RCL 03  76 *  77 STO 07  78 LBL 03  79 RCL 15  80 RCL 15  81 1  82 ST+ Z  83 -  84 RCL IND Y  85 STO 12  86 RCL IND Y  87 ST+ 12  88 -  89 STO 10  90 RCL 15

91 RCL 15  92 1  93 RCL 04  94 +  95 ST+ Z  96 -  97 RCL IND Y  98 STO 14  99 RCL IND Y 100 ST+ 14 101 - 102 STO 11 103 RCL IND 15  104 STO 09 105 ST+ X 106 ST- 12 107 ST- 14 108 RCL 03 109 X^2 110 ST/ 12 111 RCL 05 112 X^2 113 ST/ 14 114 RCL 15 115 RCL 15 116 RCL 04 117 2 118 + 119 ST+ Z 120 - 121 RCL 15 122 RCL 04 123 + 124 RCL 15 125 ST- L 126 RCL IND T 127 RCL IND T 128 RCL IND T 129 RCL IND L 130 + 131 - 132 + 133 RCL 03 134 ST+ X 135 ST/ 10

136 RCL 05 137 ST+ X 138 ST/ 11 139 * 140 / 141 STO 13 142 XEQ "dT/dX"  143 RCL 01 144 * 145 RCL IND 15 146 FS? 00 147 GTO 03 148 RDN 149 STO 09 150 RCL 15 151 RCL 15 152 RCL 04 153 1 154 + 155 ST+ Z 156 - 157 RCL 15 158 ISG X 159 CLX 160 RCL 15 161 DSE X 162 RCL IND T 163 RCL IND T 164 RCL IND T 165 RCL IND T 166 + 167 + 168 + 169 4 170 / 171 RCL 09 172 LBL 03 173 + 174 RCL 15 175 RCL 16 176 + 177 X<>Y 178 STO IND Y 179 RCL 03 180 ST- 07

181 DSE 15 182 DSE 18 183 GTO 03 184 RCL 05 185 ST- 08 186 DSE 15 187 DSE 15 188 DSE 19 189 GTO 02 190 RCL 01 191 ST+ 00 192 RCL 04 193 ENTER 194 ST+ Y 195 RCL 06           196 3 197 - 198 ST* Y 199 + 200 + 201  E3 202 / 203 RCL 04 204 + 205 23 206 + 207  E3 208 / 209 RCL 20 210 + 211 21 212 + 213 REGMOVE 214 RCL 04 215 21.02 216 + 217 STO 15 218 RCL 03 219 RCL 04 220 * 221 STO 07 222 LBL 04 223 RCL 07 224 RCL 00 225 XEQ "TXY0"

226 STO IND 15 227 RCL 03 228 ST- 07 229 DSE 15 230 GTO 04 231 RCL 20 232 20.02 233 + 234 RCL 04 235 1 236 + 237  E5 238 / 239 + 240 STO 14 241 RCL 04 242 - 243 STO 13 244 STO 15 245 RCL 05 246 RCL 06 247 * 248 STO 08 249 LBL 05 250 RCL 08 251 RCL 00 252 XEQ "TX0Z"  253 STO IND 15 254 RCL 05 255 ST- 08 256 DSE 15 257 GTO 05 258 RCL 03 259 RCL 04 260 * 261 STO 07 262 RCL 05 263 RCL 06 264 * 265 STO 08 266 LBL 06 267 RCL 08 268 RCL 00 269 XEQ "TX1Z" 270 STO IND 14

271 RCL 05 272 ST- 08 273 DSE 14 274 GTO 06 275 RCL 20 276 RCL 13 277 INT 278 DSE X 279  E3 280 / 281 + 282 20 283 + 284 STO 15 285 LBL 07 286 RCL 07 287 RCL 00 288 XEQ "TXY1"  289 STO IND 15 290 RCL 03 291 ST- 07 292 DSE 15 293 GTO 07 294 VIEW 00 295 DSE 17 296 GTO 01 297 21.02 298 RCL 20 299  E3 300 / 301 + 302 RCL 04 303 1 304 + 305  E5 306 / 307 + 308 RCL 00 309 RTN 310 GTO 10 311 END

 
       ( 473 bytes / 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 - y² ¶²T/¶y² + 2 y.z ¶²T/¶y¶z - z² ¶²T/¶z²

  such that  x₀ = 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
 

6°) 2D Poisson Equation
 

 "2DPOIS" solves     A ¶²U/¶x² + B ¶²UT/¶x¶y + C ¶²U/¶y² = 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"
 
 

01 LBL "2DPOIS"  02 RCL 02  03 1  04 +  05 STO 00  06  E3  07 STO 10  08 /  09 16.015  10 STO 11  11 +  12 STO 09  13 CLX  14 STO 05  15 LBL 01  16 RCL 05  17 XEQ "UX0"  18 STO IND 09  19 RCL 01  20 ST+ 05  21 ISG 09  22 GTO 01  23 RCL 00  24 1  25 %  26 RCL 00  27 RCL 04  28 LASTX  29 +  30 *  31 +  32 RCL 10  33 /  34 RCL 11  35 +  36 STO 09  37 CLX  38 STO 06

39 LBL 02  40 RCL 06  41 XEQ "U0Y"  42 STO IND 09     43 RCL 03  44 ST+ 06  45 ISG 09  46 GTO 02  47 DSE 09  48 CLX  49 RCL 09  50 STO 08  51 RCL 10  52 *  53 INT  54 RCL 10  55 /  56 STO 09  57 CLX  58 STO 05  59 LBL 03  60 RCL 05  61 XEQ "UX1"  62 STO IND 09  63 RCL 01  64 ST+ 05  65 ISG 09  66 GTO 03  67 RCL 08  68 FRC  69 RCL 02  70 +  71 16  72 +  73 STO 09  74 CLX  75 STO 06  76 LBL 04

77 RCL 06  78 XEQ "U1Y"  79 STO IND 09     80 RCL 03  81 ST+ 06  82 ISG 09  83 GTO 04  84 LBL 05  85 CLX  86 STO 15  87 RCL 02  88 DSE X  89 STO 07  90 RCL 01  91 *  92 STO 05  93 RCL 04  94 DSE X  95 STO 08  96 RCL 03  97 *  98 STO 06  99 RCL 00 100 RCL 04 101 * 102 14 103 + 104 STO 09 105 LBL 06 106 RCL 06 107 RCL 05 108 XEQ "AXY" 109 RCL 01 110 X^2 111 / 112 STO 10 113 RCL 06 114 RCL 05

115 XEQ "BXY" 116 RCL 01 117 RCL 03 118 * 119 ST+ X 120 / 121 STO 11 122 RCL 06 123 RCL 05 124 XEQ "CXY"   125 RCL 03 126 X^2 127 / 128 STO 12 129 RCL 06 130 RCL 05 131 XEQ "DXY" 132 RCL 01 133 ST+ X 134 / 135 STO 13 136 RCL 06 137 RCL 05 138 XEQ "EXY" 139 RCL 03 140 ST+ X 141 / 142 STO 14 143 RCL 09 144 1 145 + 146 RCL IND X 147 RCL 10 148 RCL 13 149 - 150 RCL 11 151 ST- 13 152 -

153 * 154 RCL 09 155 DSE X 156 X<>Y 157 RCL IND Y 158 RCL 10 159 RCL 13 160 + 161 * 162 + 163 STO 13 164 RCL 09 165 RCL 00 166 - 167 RCL IND X   168 RCL 12 169 RCL 14 170 + 171 RCL 11 172 ST+ 14 173 - 174 * 175 ST+ 13 176 RCL 09 177 RCL 00 178 + 179 RCL IND X 180 RCL 12 181 RCL 14 182 - 183 * 184 ST+ 13 185 RCL 06 186 RCL 05 187 XEQ "GXY" 188 ST- 13 189 RCL 09 190 RCL 00

191 1 192 + 193 + 194 RCL 09 195 LASTX 196 - 197 RCL IND Y 198 RCL IND Y 199 + 200 RCL 11 201 * 202 ST+ 13 203 RCL 06 204 RCL 05 205 XEQ "FXY"   206 RCL 10 207 RCL 11 208 - 209 RCL 12 210 + 211 ST+ X 212 + 213 RCL 13 214 X<>Y 215 / 216 RCL IND 09 217 - 218 PI 219 LN 220 FS? 00 221 SIGN 222 * 223 ST+ IND 09 224 ABS 225 RCL 15 226 X<Y? 227 X<>Y 228 STO 15

229 RCL 01 230 ST- 05 231 DSE 09 232 DSE 07 233 GTO 06 234 DSE 09 235 DSE 09 236 RCL 03 237 ST- 06 238 RCL 02 239 DSE X 240 STO 07 241 RCL 01 242 * 243 STO 05           244 DSE 08 245 GTO 06 246 RCL 15 247 VIEW X 248 RND 249 X#0? 250 GTO 05 251 RCL 00 252 1 253 % 254 RCL 00 255 RCL 04 256 LASTX 257 + 258 * 259 + 260  E3 261 / 262 16.015 263 + 264 CLD 265 END

 
      ( 388 bytes / SIZE 016+(M+1).(N+1) )
 
 

STACK	INPUTS	OUTPUTS
X	/	16.eee.N+1

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

Example:       ¶²U/¶x² + ¶²UT/¶x¶y + ¶²U/¶y² = -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
 

7°)  Wave Equation ( Explicit Method )
 

-"WAVE"  solves equations of the type:     ¶²W/¶x² = A(x,y) ¶²W/¶y² + 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 = x₀                                    ( 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"
 
 

01 LBL "WAVE"  02 RCL 04  03 ENTER  04 STO 10  05 14  06 +  07 STO 09  08 +  09 STO 11  10 ISG 11  11 LBL 00  12 RCL 03  13 RCL 10  14 *  15 STO 05  16 XEQ "W0Y"  17 STO IND 09  18 RCL 05  19 XEQ "dW0Y"  20 RCL 01  21 *  22 RCL IND 09  23 +  24 STO IND 11  25 DSE 09  26 DSE 11  27 DSE 10  28 GTO 00  29 CLX

30 XEQ "W0Y"   31 STO 14  32 1  33 STO 08  34 RCL 02  35 STO 06  36 GTO 14  37 LBL 01  38 RCL 02  39 STO 06  40 LBL 02  41 RCL 04  42 ST+ X  43 15  44 +  45 RCL 08  46 -  47 .1  48 %  49 +  50 RCL 04  51 -  52 ISG X  53 RCL 08  54 +  55 STO 07  56 INT  57 RCL 04  58 15

59 +  60 -  61 RCL 03  62 *  63 STO 05  64 LBL 03  65 RCL 05  66 RCL 00  67 XEQ "AXY"   68 RCL 03  69 X^2  70 /  71 STO 09  72 RCL 05  73 RCL 00  74 XEQ "DXY"  75 RCL 09  76 RCL 01  77 X^2  78 1/X  79 -  80 ST+ X  81 -  82 RCL IND 07  83 *  84 STO 10  85 RCL 05  86 RCL 00  87 XEQ "CXY"

88 RCL 03  89 ST+ X  90 /  91 RCL 09  92 X<>Y  93 ST- 09  94 +  95 RCL 07  96 1  97 +  98 X<>Y  99 RCL IND Y  100 * 101 ST+ 10 102 RCL 07 103 1 104 - 105 RCL IND X   106 RCL 09 107 * 108 ST+ 10 109 RCL 05 110 RCL 00 111 XEQ "BXY" 112 RCL 01 113 ST+ X 114 / 115 CHS 116 STO 09

117 RCL 01 118 X^2 119 1/X 120 ST+ 09 121 - 122 RCL 07 123 RCL 04 124 - 125 1 126 - 127 X<>Y 128 RCL IND Y 129 * 130 ST+ 10 131 RCL 05 132 RCL 00 133 XEQ "EXY"  134 ST+ 10 135 RCL 04 136 RCL 07 137 + 138 1 139 + 140 RCL 10 141 RCL 09 142 / 143 STO IND Y 144 RCL 03 145 ST+ 05

146 ISG 07 147 GTO 03 148 RCL 04 149 ST+ X 150  E3 151 / 152 14 153 + 154  E3 155 / 156 15 157 + 158 RCL 04 159 + 160 REGMOVE  161 LBL 14 162 RCL 01 163 ST+ 00 164 1 165 ST+ 08 166 VIEW 00 167 DSE 06 168 GTO 02 169 RCL 04 170 8 171 + 172 ST+ X 173 RCL 08 174 -

175  E3 176 / 177 RCL 04 178 RCL 08 179 + 180 14 181 + 182 + 183 STO Y 184 RCL 04 185 .1 186 % 187 + 188 - 189 2 190 - 191 X<>Y 192 RCL 08         193 1 194 - 195 RCL 03 196 * 197 RCL 00 198 RTN 199 GTO 01 200 END

 
     ( 295 bytes / 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:     ¶²W/¶x² =  ¶²W/¶y² - y ¶W/¶x + y ¶W/¶y - x² W + x.y exp(-x.y)

    with  W(0,y) = 1  and  ¶W/¶x(0,y) = - y      (  x₀ = 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 x₀ = 0   h = 0.05  M =10   k = 0.05  N = 20   you'll get  W(0.5,0.5) ~ 0.7750
-With x₀ = 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
 

8°)  2D Wave Equation ( Explicit Method )
 

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

    ¶²W/¶x² = A(x,y,z) ¶²W/¶y² + B(x,y,z) ¶²W/¶y¶z + C(x,y,z) ¶²W/¶z² + 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 = x₀                                    ( 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"
 
 

01 LBL "2DWAVE"  02 RCL 04  03 1  04 +  05 STO 09  06 RCL 06  07 LASTX  08 +  09 STO 10  10 *  11 STO 14  12 20  13 +  14 STO 13  15 RCL 03  16 RCL 04  17 *  18 STO 07  19 RCL 05  20 RCL 06  21 *  22 STO 08  23 LBL 00  24 RCL 08  25 RCL 07  26 XEQ "W0YZ"  27 STO IND 13  28 RCL 03  29 ST- 07  30 DSE 13  31 DSE 09  32 GTO 00  33 RCL 03  34 RCL 04  35 *  36 STO 07  37 RCL 05  38 ST- 08  39 RCL 04  40 1  41 +  42 STO 09  43 DSE 10  44 GTO 00  45 RCL 14  46 RCL 04  47 -  48 18  49 +  50 STO 17  51 RCL 14  52 RCL 04  53 DSE X  54 RCL 06  55 DSE X  56 *  57 +  58 20  59 +  60 STO 18  61 RCL 06

62 DSE X  63 STO 10  64 RCL 05  65 *  66 STO 08  67 LBL 01  68 RCL 04  69 DSE X  70 STO 09  71 RCL 03  72 *  73 STO 07  74 LBL 02  75 RCL 08  76 RCL 07  77 XEQ "dW0YZ"  78 RCL 01  79 *  80 RCL IND 17  81 +  82 STO IND 18  83 RCL 03  84 ST- 07  85 DSE 18  86 DSE 17  87 DSE 09  88 GTO 02  89 DSE 17  90 DSE 17  91 RCL 05  92 ST- 08  93 DSE 10  94 GTO 01  95 1  96 STO 14  97 RCL 02  98 STO 20  99 GTO 14 100 LBL 13 101 RCL 02 102 STO 20 103 LBL 10 104 RCL 06 105 RCL 14 106 ST+ X 107 - 108 STO 12 109 5 110 + 111 RCL 04 112 RCL 14 113 ST+ X 114 - 115 STO 11 116 5 117 + 118 * 119 20 120 + 121 STO 17 122 STO 18

123 STO 19 124 RCL 11 125 5 126 + 127 ST+ X 128 2 129 + 130 ST- 17 131 RCL 12 132 3 133 + 134 RCL 11 135 LASTX 136 + 137 * 138 STO 15            139 ST+ 18 140 ST+ 19 141 LASTX 142 1 143 + 144 ST- 18 145 RCL 11 146 LASTX 147 + 148 RCL 12 149 LASTX 150 + 151 * 152 ST+ 15 153 ST+ 19 154 RCL 06 155 RCL 14 156 - 157 STO 10 158 RCL 05 159 * 160 STO 08 161 RCL 14 162 ST- 10 163 ISG 10 164 LBL 11 165 RCL 04 166 RCL 14 167 - 168 STO 09 169 RCL 03 170 * 171 STO 07 172 RCL 14 173 ST- 09 174 ISG 09 175 LBL 12 176 RCL 08 177 RCL 07 178 RCL 00 179 XEQ "AXYZ" 180 RCL 03 181 X^2 182 / 183 STO 11

184 RCL 08 185 RCL 07 186 RCL 00 187 XEQ "CXYZ" 188 RCL 05 189 X^2 190 / 191 STO 12 192 RCL 08 193 RCL 07 194 RCL 00 195 XEQ "GXYZ"  196 RCL 01 197 X^2 198 1/X 199 RCL 12 200 - 201 RCL 11 202 - 203 ST+ X 204 + 205 RCL IND 18 206 * 207 STO IND 19 208 RCL 08 209 RCL 07 210 RCL 00 211 XEQ "HXYZ" 212 ST+ IND 19 213 RCL 08 214 RCL 07 215 RCL 00 216 XEQ "EXYZ" 217 RCL 03 218 ST+ X 219 / 220 STO 13 221 RCL 18 222 1 223 + 224 RCL IND X 225 RCL 11 226 RCL 13 227 + 228 * 229 ST+ IND 19 230 RCL 18 231 DSE X 232 RCL IND X 233 RCL 11 234 RCL 13 235 - 236 * 237 ST+ IND 19 238 RCL 08 239 RCL 07 240 RCL 00 241 XEQ "FXYZ" 242 RCL 05 243 ST+ X 244 /

245 STO 13 246 RCL 18 247 RCL 04 248 + 249 RCL 14 250 ST+ X 251 - 252 3 253 + 254 RCL IND X     255 RCL 12 256 RCL 13 257 + 258 * 259 ST+ IND 19 260 RCL 18 261 RCL 04 262 - 263 RCL 14 264 ST+ X 265 + 266 3 267 - 268 RCL IND X 269 RCL 12 270 RCL 13 271 - 272 * 273 ST+ IND 19 274 RCL 08 275 RCL 07 276 RCL 00 277 XEQ "BXYZ" 278 RCL 03 279 RCL 05 280 * 281 4 282 * 283 / 284 STO 12 285 RCL 18 286 RCL 04 287 4 288 + 289 RCL 14 290 ST+ X 291 - 292 + 293 RCL 18 294 LASTX 295 - 296 RCL 04 297 RCL 14 298 DSE X 299 ST+ X 300 - 301 STO L 302 RCL 18 303 ST+ Y 304 ST- L 305 RCL IND T

306 RCL IND T 307 RCL IND T 308 RCL IND L 309 + 310 - 311 + 312 RCL 12 313 * 314 ST+ IND 19 315 RCL 08 316 RCL 07 317 RCL 00 318 XEQ "DXYZ"  319 RCL 01 320 ST+ X 321 / 322 STO 12 323 RCL 01 324 X^2 325 1/X 326 + 327 RCL IND 17 328 * 329 ST- IND 19 330 RCL 01 331 X^2 332 1/X 333 RCL 12 334 - 335 ST/ IND 19 336 DSE 19 337 DSE 18 338 DSE 17 339 DSE 09 340 GTO 12 341 RCL 05 342 ST- 08 343 2 344 ST- 18 345 X^2 346 ST- 17 347 DSE 10 348 GTO 11 349 RCL 04 350 RCL 14 351 ST+ X 352 5 353 - 354 - 355 RCL 06 356 LASTX 357 - 358 * 359 21.021 360 + 361 RCL 15 362  E6 363 / 364 + 365 REGMOVE 366 LBL 14

367 1 368 ST+ 14 369 RCL 01 370 ST+ 00 371 VIEW 00 372 DSE 20 373 GTO 10 374 RCL 04 375 RCL 14 376 ST+ X 377 5 378 - 379 STO 11            380 - 381 STO 15 382 RCL 06 383 LASTX 384 - 385 * 386 STO 12 387  E3 388 / 389 21.02 390 + 391 RCL 15 392  E5 393 / 394 + 395 STO 13 396 2 E-5 397 - 398 RCL 12 399 + 400 RCL 04 401 RCL 11 402 2 403 + 404 - 405 RCL 06 406 LASTX 407 - 408 * 409  E3 410 / 411 + 412 STO 12 413 RCL 05 414 RCL 03 415 RCL 14 416 DSE X 417 ST* Z 418 * 419 RCL 00 420 RCL 13 421 SIGN 422 CLX 423 RCL 12 424 RDN 425 RTN 426 GTO 13 427 END

 
     ( 619 bytes / 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:     ¶²W/¶x² = ¶²W/¶y² + 2 ¶²W/¶y¶z + ¶²W/¶z² - ¶W/¶x - y ¶W/¶y + z ¶W/¶z - ( y + z )² W + 3 exp(-x-y.z)

    x₀ = 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
 

9°)  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 = §_a^b  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 = §_a^b  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.
 
 
 

01 LBL "SHOOT"  02 STO 06  03 X<>Y  04 STO 05  05 RCL 00  06 2  07 MOD  08 ST+ 00  09 RCL 02  10 RCL 01  11 -  12 RCL 00  13 ST+ X  14 /  15 STO 10  16 SF 10  17 RCL 06  18 XEQ 02  19 STO 14  20 LBL 01  21 VIEW 05  22 RCL 05  23 XEQ 02  24 ENTER  25 ENTER  26 X<> 14  27 -  28 X#0?  29 /  30 RCL 06

31 RCL 05  32 STO 06  33 STO T  34 -  35 *  36 +  37 STO 05            38 -  39 RND  40 X#0?  41 GTO 01  42 2  43 STO 15  44 X^2  45 STO 18  46 20  47 STO 19  48 RCL 05  49 CF 10  50 LBL 02  51 STO 09  52 RCL 01  53 STO 07  54 RCL 00  55 STO 16  56 RCL 03  57 STO 08  58 FS? 10  59 GTO 03  60 RCL 09

61 X<>Y  62 STO 20  63 RCL 07  64 XEQ "L"  65 STO 17  66 LBL 03  67 RCL 09            68 RCL 08  69 RCL 07  70 XEQ "T"  71 RCL 10  72 ST+ 07  73 *  74 STO 12  75 RCL 09  76 +  77 STO 13  78 LASTX  79 RCL 10  80 *  81 STO 11  82 RCL 08  83 +  84 RCL 07  85 XEQ "T"  86 RCL 10  87 ST* 13  88 *  89 ST+ 12  90 ST+ 12

91 RCL 09  92 +  93 ENTER  94 X<> 13  95 ST+ 11  96 ST+ 11  97 RCL 08  98 +  99 RCL 07 100 XEQ "T" 101 RCL 10          102 ST+ 07 103 ST+ X 104 ST* 13 105 * 106 ST+ 12 107 RCL 09 108 + 109 ENTER 110 X<> 13 111 ST+ 11 112 RCL 08 113 + 114 RCL 07 115 XEQ "T" 116 RCL 10 117 ST* 13 118 * 119 RCL 12 120 +

121 3 122 / 123 ST+ 09 124 RCL 13 125 RCL 11 126 + 127 3 128 / 129 ST+ 08 130 FS? 10 131 GTO 03 132 ISG 19 133 CLX 134 RCL 09 135 RCL 08 136 FC? 00 137 STO IND 19  138 RCL 07 139 XEQ "L" 140 RCL 15 141 X<> 18 142 STO 15 143 * 144 ST+ 17 145 LBL 03 146 DSE 16 147 GTO 03 148 RCL 08 149 RCL 04 150 -

151 FS? 10 152 RTN 153 STO 14 154 CLX 155 2 156 / 157 ST- 17 158 RCL 19 159  E3 160 / 161 20 162 + 163 FS? 00 164 CLX 165 STO 08 166 RCL 17           167 RCL 10 168 * 169 1.5 170 / 171 STO 07 172 RCL 14 173 SIGN 174 CLX 175 RCL 09 176 STO 06 177 RCL 05 178 END

 
     ( 255 bytes / SIZE 021+N or 021 )
 
 

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⁺¹  y ( 1 + y' ² ) ^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' ² ) / 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⁺¹  y ( 1 + y' ² ) ^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   §_a^b  T( x , y , y' ) = y'(b) - y'(a)
 

10°)  Euler-Lagrange Equation
 

"EULAG" finds the minimum - or extremum - of the integral   I = §_a^b  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 = ¶²L/¶x¶y' + y'  ¶²L/¶y¶y' + y'' ¶²L/¶y'²

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

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

     y(x+h) ~ y(x) + h.y'(x) + ( h² / 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.
 
 

01 LBL "EULAG"  02 STO 06  03 X<>Y  04 STO 05  05 RCL 00  06 2  07 MOD  08 ST+ 00  09 RCL 02  10 RCL 01  11 -  12 RCL 00  13 /  14 STO 10  15 SF 10  16 RCL 06  17 XEQ 02  18 STO 11  19 LBL 01  20 VIEW 05  21 RCL 05  22 XEQ 02  23 ENTER  24 ENTER  25 X<> 11  26 -  27 X#0?  28 /  29 RCL 06  30 RCL 05  31 STO 06  32 STO T  33 -  34 *  35 +  36 STO 05

37 -  38 RND  39 X#0?  40 GTO 01  41 20  42 STO 19  43 2  44 STO 16            45 X^2  46 STO 15  47 RCL 05  48 CF 10  49 LBL 02  50 STO 09  51 RCL 01  52 STO 07  53 RCL 00  54 STO 18  55 RCL 03  56 STO 08  57 FS? 10  58 GTO 03  59 RCL 09  60 X<>Y  61 STO 20  62 RCL 07  63 XEQ "L"  64 CHS  65 STO 14  66 LBL 03  67 RCL 09  68 RCL 08  69 RCL 10  70 +  71 RCL 07  72 XEQ "L"

73 STO 12  74 RCL 09  75 RCL 08  76 RCL 10  77 -  78 RCL 07  79 XEQ "L"  80 ST- 12  81 RCL 10            82 ST+ X  83 ST/ 12  84 RCL 09  85 RCL 08  86 RCL 07  87 RCL 10  88 ST+ T  89 +  90 XEQ "L"  91 STO 13  92 RCL 09  93 RCL 08  94 RCL 07  95 RCL 10  96 ST- T  97 -  98 XEQ "L"  99 ST+ 13 100 RCL 09 101 RCL 08 102 RCL 07 103 RCL 10 104 ST+ T 105 - 106 XEQ "L" 107 ST- 13 108 RCL 09

109 RCL 08 110 RCL 07 111 RCL 10 112 ST- T 113 + 114 XEQ "L" 115 ST- 13 116 RCL 13 117 RCL 10 118 ST+ X 119 X^2 120 STO 17 121 / 122 ST- 12 123 RCL 09 124 RCL 08 125 RCL 10 126 ST+ Z 127 + 128 RCL 07 129 XEQ "L" 130 STO 13 131 RCL 09 132 RCL 08 133 RCL 10 134 ST- Z 135 - 136 RCL 07 137 XEQ "L" 138 ST+ 13 139 RCL 09 140 RCL 08 141 RCL 10 142 ST+ Z 143 - 144 RCL 07

145 XEQ "L" 146 ST- 13 147 RCL 09 148 RCL 08 149 RCL 10 150 ST- Z 151 + 152 RCL 07 153 XEQ "L" 154 ST- 13 155 RCL 13          156 RCL 17 157 / 158 RCL 09 159 * 160 ST- 12 161 RCL 09 162 RCL 10 163 + 164 RCL 08 165 RCL 07 166 XEQ "L" 167 STO 13 168 RCL 09 169 RCL 10 170 - 171 RCL 08 172 RCL 07 173 XEQ "L" 174 ST+ 13 175 RCL 09 176 RCL 08 177 RCL 07 178 XEQ "L" 179 STO 17 180 ST+ X

181 ST- 13 182 RCL 13 183 RCL 10 184 X^2 185 / 186 ST/ 12 187 RCL 12 188 RCL 10 189 * 190 2 191 / 192 RCL 09          193 + 194 RCL 10 195 * 196 ST+ 08 197 RCL 12 198 RCL 10 199 ST+ 07 200 * 201 ST+ 09 202 FS? 10 203 GTO 03 204 RCL 17 205 RCL 15 206 X<> 16 207 STO 15 208 * 209 ST+ 14 210 ISG 19 211 CLX 212 RCL 08 213 FC? 00 214 STO IND 19 215 LBL 03 216 DSE 18

217 GTO 03 218 RCL 08 219 RCL 04 220 - 221 FS? 10 222 RTN 223 STO 11 224 RCL 09 225 RCL 08          226 RCL 07 227 XEQ "L" 228 ST+ 14 229 RCL 19 230  E3 231 / 232 20 233 + 234 FS? 00 235 CLX 236 STO 08 237 RCL 14 238 RCL 10 239 * 240 3 241 / 242 STO 07 243 RCL 11 244 SIGN 245 CLX 246 RCL 09 247 STO 06 248 RCL 05 249 END

 
        ( 354 bytes / SIZE 021+N or 021 )
 
 

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⁺¹  y ( 1 + y' ² ) ^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⁺¹  y ( 1 + y' ² ) ^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