HP35S - Ordinary Differential Equations (and incidentally Lambert W Function)

[Roberto Volpi]

Hi all

The following program, using the Runge-Kutta algorithm, can compute numerical values of an ordinary differential equation (ODE), both 1st and 2nd order, given initial conditions.

LINE INSTRUCTION
O001 LBL O
O002 INPUT X
O003 STO E
O004 INPUT Y
O005 STO W
O006 F2?
O007 XEQ O083
O008 INPUT U
O009 INPUT H
O010 2
O011 STO F
O012 3
O013 STO G
O014 RCL E
O015 STO T
O016 RCL W
O017 STO V
O018 F2?
O019 XEQ O086
O020 XEQ F001
O021 RCL xH
O022 F2?
O023 RCL /F
O024 STO A
O025 RCL /F
O026 F2?
O027 XEQ O089
O028 RCL +W
O029 STO V
O030 RCL H
O031 RCL /F
O032 RCL +E
O033 STO T
O034 F2?
O035 XEQ O093
O036 XEQ F001
O037 RCL xH
O038 F2?
O039 RCL /F
O040 STO B
O041 F2?
O042 GTO O097
O043 RCL/F
O044 RCL +W
O045 STO V
O046 XEQ F001
O047 RCL xH
O048 F2?
O049 RCL /F
O050 STO C
O051 F2?
O052 XEQ O100
O053 RCL +W
O054 STO V
O055 RCL E
O056 RCL +H
O057 STO T
O058 XEQ F001
O059 RCL xH
O060 F2?
O061 RCL /F
O062 STO D
O063 RCL H
O064 STO +E
O065 F2?
O066 XEQ O107
O067 RCL D
O068 RCL +A
O069 RCL +B
O070 RCL +B
O071 RCL +C
O072 RCL +C
O073 RCL /G
O074 RCL /F
O075 F2?
O076 GTO O115
O077 STO +W
O078 RCL E
O079 RCL U
O080 X=Y?
O081 GTO O118
O082 GTO O014
O083 INPUT Z
O084 STO Q
O085 RTN
O086 RCL Q
O087 STO S
O088 RTN
O089 RCL +Q
O090 RCL xH
O091 RCL /F
O092 RTN
O093 RCL A
O094 RCL +Q
O095 STO S
O096 RTN
O097 RCL +Q
O098 STO S
O099 GTO O046
O100 RCL xF
O101 RCL +Q
O102 STO S
O103 RCL C
O104 RCL +Q
O105 RCL xH
O106 RTN
O107 RCL A
O108 RCL +B
O109 RCL +C
O110 RCL /G
O111 RCL +Q
O112 RCL xH
O113 STO +W
O114 RTN
O115 lastX
O116 STO +Q
O117 GTO O078
O118 VIEW W
O119 VIEW Q
O120 GTO O118


It is possible to choose the magnitude of h: the smaller h, the better approximation will be reached, but at cost of a longer time requested by the calculator.

The ODE to be solved will be stored in LBL F in RPN form and using the following registers:
T=x; V=y; S=y’

1st Example: y’=x+y
LBL F
RCL T
RCL +V
RTN

2nd Example: y’’=xy’+y will be:
LBL F
RCL T
RCL xS
RCL +V
RTN

As you can see, only right-hand side of the equation is entered. Our HP35S will consider it a 1st order ODE by default, and a 2nd order ODE after setting FLAG 2 before executing the program.

Remark: as, to speed calculations, register F contains already 2 and register G contains already 3, you are strongly encouraged to take advantage of it by entering the ODE in LBL F in case the value 2, or 3, are part of the ODE. As an example, the ODE y’=2xy will be entered as:

LBL F
RCL T
RCL xV
RCL xF
RTN

Instructions to solve a 1st order ODE already in LBL F:

XEQ O
> INPUT X: digit initial value of X (x0); press R/S
> INPUT Y: digit initial value of Y (y0); press R/S
> INPUT U: digit value of X (x1) for which we would find correspondent y1 value; press R/S
> INPUT H: digit value of h. A good approximation is already possible with 0.1. Then press R/S
> W=y1; the calculator will display y1 as value of W register.

Instructions to solve a 2st order ODE:

Same as above, but:

- SET FLAG 2, so the function stored in LBL F will be read as function of y’’ instead of function of y’.
- After INPUT U, the calculator will prompt INPUT Z, so we digit y’0 value, and then R/S
- The calculator will give w=y1, press R/S, and we obtain Q=y’1 value. By pressing R/S again you can see again the value stored in W, press again R/S to see Q, and again ad libitum.


USE THIS PROGRAM TO APPROXIMATE THE LAMBERT W FUNCTION

Lambert W function is defined as inverse of y=xe^x

In implicit form, it can be expressed as x=ye^y, but it is not very useful to use it this way.

We can obtain a decent approximation with this program by making use of ODE of Lambert W function y’=1/e^y(1+y), so we can input it as follows:

LBL F
RCL V
e^x
RCL xV
lastX
+
1/x
RTN

Initial values x0=0 and y0=0, so now we try to compute some values and estimate how long will our HP35S take, and finally compare results with those given by wolframalpha.com

Lambert W of x=1

W=0.578265978405, after 1 second and setting h=1
W=0.5671439206, after 7 seconds and setting h=0.1
W=0.567143290409… with wolframalpha.com

Lambert W of x=7

W=1.5270617364, after 5 second and setting h=1
W=1.52434535808, after 47 seconds and setting h=0.1
W=1.52434520848… with wolframalpha.com

I hope no typo has been made…

Enjoy!