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/¶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
| 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. |
-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/¶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
¶2T/¶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
¶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/¶x¶y + 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/¶y¶z + C ¶2T/¶z2 + D ¶2T/¶x¶y + E ¶2T/¶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)
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/¶y¶z + ¶2T/¶z2 + y ¶2T/¶x¶y - z ¶2T/¶x¶z - 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/¶y¶z , ¶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/¶y¶z
• 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/¶y¶z
, ¶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/¶y¶z - 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
"2DPOIS" solves A ¶2U/¶x2 + B ¶2UT/¶x¶y + 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/¶x¶y + ¶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/¶x2 = ¶2W/¶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/¶y¶z + 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/¶y¶z + ¶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
"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)
"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/¶x¶y' + y' ¶2L/¶y¶y' + 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