HP35S  Ordinary Differential Equations (and incidentally Lambert W Function)

01182024, 01:21 PM
(This post was last modified: 01182024 01:27 PM by Roberto Volpi.)
Post: #1




HP35S  Ordinary Differential Equations (and incidentally Lambert W Function)
Hi all
The following program, using the RungeKutta 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 righthand 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! Put a calculator into your life! 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)