(71B) EulerTaylor method for the HP71B

12072019, 04:36 PM
(This post was last modified: 12112019 10:04 PM by Namir.)
Post: #1




(71B) EulerTaylor method for the HP71B
EulerTaylor method
=================== To integrate f'(x,y) from x=a to x=b, give y(a), we use an extended version of the Euler method that includes the first AND second derivatives. Since we know f'(x,y) we can derive f''(x,y) and use the following formula: y = y + h*f'(x,y) + h^2/2*f''(x,y) HP71B Listing ============== Code: 10 REM EULER METHOD WITH 2 DERIVATIVES Usage ===== 1) With the values in the DATA statements being the ones you want and the definitions of functions FNE, FNX, and FND representing your problem to solve, just press the [RUN] button. The program displays "WAIT PLEASE ..." as it calculates y at x=b. 2) Displays the value of y at x=b. Press [f][Cont] to resume. 3) Displays the exact value of y at x=b. Press [f][Cont] to resume. 4) Displays the % error between the calculated and exact values of y at x=b. Code Customization ================== 1) Edit the DATA statement to use different values for a, b, y(a), and h. 2) Edit the definitions of functions FNE, FNX, and FND to use a different exact function, a different first derivative, and a different second derivatve, respectively. In the case of FNE, edit the DEF FNE(X) statement to implement the eaxct solutions as a function of variable x only. If that solution is not available use something like: DEF FNE(X)=1E99 HP71B Listing 2 ================ This version should be faster since the FOR loop calculates the values of the first and second derivative directly and does not resort to calling functions defined in DEF statements. Code: 10 REM EULER METHOD WITH 2 DERIVATIVES ===== 1) With the values in the DATA statements being the ones you want and the definitions of functions FNE, FNX, and FND representing your problem to solve, just press the [RUN] button. The program displays "WAIT PLEASE ..." as it calculates y at x=b. 2) Displays the value of y at x=b. Press [f][Cont] to resume. 3) Displays the exact value of y at x=b. Press [f][Cont] to resume. 4) Displays the % error between the calculated and exact values of y at x=b. Code Customization ================== 1) Edit the DATA statement to use different values for a, b, y(a), and h. 2) Edit line 90 to change how the first derivative is calculated and assigned to variable D1. 3) Edit line 100 to change how the first derivative is calculated and assigned to variable D2. You can use variable D1 when the value of teh first derivative is needed. 4) Edit the DEF FNE(X) statement to implement the eaxct solutions as a function of variable x only. If that solution is not available use something like: DEF FNE(X)=1E99 

12102019, 06:02 PM
Post: #2




RE: (71B) EulerTaylor method for the HP71B
Hi Namir, just two comments:
1.) at relative error calculation the divisor is Y0 not 10, I guess 2.) I am not familiar with 71B, but in the function definitions there is not important the case sensitivy (y or Y)? +1) you know that, there is a cheating when you defined the second derivative, I mean a differentiation procedure is required, because thats why we use computers... Csaba 

12112019, 10:07 PM
Post: #3




RE: (71B) EulerTaylor method for the HP71B
(12102019 06:02 PM)Csaba Tizedes Wrote: Hi Namir, just two comments: 1) Thatls for spotting the error. I corrected it in both listings. 2) HP71B BASIC, as mots BASICs, is not case sensitive. 3) It would be nice if someone can numerically approximate the d2y/dx2 for dy/dx = f(x,y). Any idea? When say dy/dx=y, approximating the second derivative with x and an increment of x leads no where! I am open to suggestions. Namir 

12132019, 02:12 PM
(This post was last modified: 12142019 07:11 PM by Albert Chan.)
Post: #4




RE: (71B) EulerTaylor method for the HP71B
(12112019 10:07 PM)Namir Wrote: 3) It would be nice if someone can numerically approximate the d2y/dx2 for dy/dx = f(x,y). Any idea? y'' ≈ Δy' / Δx = Δy' / h Δy' ≈ y'(x+h, y + h y'(x,y))  y'(x,y) → y + h y' + ½ h² y'' ≈ y + h y' + ½ h Δy' 10 DEF FND(X,Y)=X*Y 20 X=0 @ Y=1 @ X1=1 @ H=.01 30 N=INT((X1X)/H+.5) 40 FOR I=1 TO N 50 D1=FND(X,Y) 60 X=X+H @ Y=Y+H*D1 70 Y=Y+H/2*(FND(X,Y)D1) 80 NEXT I 90 Y1=EXP(.5*X*X) 100 DISP "y =";Y 110 DISP "exact=";Y1 120 DISP "%err =";100*(YY1)/Y1 >RUN y = 1.64871423741 exact = 1.6487212707 %err =4.26590602365E4 Edit:To compare apples to apples, I removed the assumption final X=X1 So, both exact and estimate uses final X as endpoint. line 90: Y1=EXP(.5*X1*X1) changed to Y1=EXP(.5*X*X) 

12132019, 04:07 PM
Post: #5




RE: (71B) EulerTaylor method for the HP71B
Thank you very much Albert for your approximation and the updated HP71B listing. My hats off to you sir!!
:) Namir 

12132019, 04:10 PM
(This post was last modified: 12132019 04:11 PM by Namir.)
Post: #6




RE: (71B) EulerTaylor method for the HP71B
(12132019 02:12 PM)Albert Chan Wrote:(12112019 10:07 PM)Namir Wrote: 3) It would be nice if someone can numerically approximate the d2y/dx2 for dy/dx = f(x,y). Any idea? Do I have your permission to update/adapt my HP71B listing using your approximation for the second derivative? I will give you credit for that update (and also for the update of the HP41C Listing). Namir 

12132019, 04:43 PM
(This post was last modified: 12132019 04:44 PM by Albert Chan.)
Post: #7




RE: (71B) EulerTaylor method for the HP71B
(12132019 04:10 PM)Namir Wrote: Do I have your permission to update/adapt my HP71B listing using your approximation for the second derivative? Sure. Glad it helps. FYI, you can delete line for y 2nd order correction, and see how much y'' helps. With the same example, and *only* first derivative: %error goes up from 4.266e4 to to 0.6612 2nd order correction improved accuracy 1550X 

12132019, 09:16 PM
Post: #8




RE: (71B) EulerTaylor method for the HP71B
Albert,
For last last month I have been studying numerical methods for ODE that are based on RungeKutta methods. There is a VAST number of variants of the RK methods! I happen to browse at early pages of a book by Leon Lapidus and I noticed the Tayler expansion. At that moment, something registered in my brain and I decided to extend the basic Euler method to include the second derivative. The increase in accuracy relative to the CPU effort is very good when you compare it with different variants of the RungeKutta methods that calculate k1 through k5 or k6 intermediate values per iteration. I realized that extending the Euler method is very suitable for our beloved calculator. Your contibution to approximate the second derivative is key in making the method more practical to use. Namir 

12142019, 03:33 AM
Post: #9




RE: (71B) EulerTaylor method for the HP71B
Adapting my original HP71B to Albert Chan's approximaton of the second derivative, we get:
Code: 10 REM EULER METHOD WITH 2 DERIVATIVES 

12142019, 08:44 AM
Post: #10




RE: (71B) EulerTaylor method for the HP71B
(12132019 02:12 PM)Albert Chan Wrote: Δy' ≈ y'(x+h, y + h y'(x,y))  y'(x,y) Thanks. You're lucky man, you have more spare time than me. It is an endyearrush in power industry and today (Saturday) is a work day in Hungary (in this December two Saturdays are workdays  these are "relocated" to increase the Christmas holiday period). So, you're lucky and thanks for the substitutions and the result  is looks like a PredictorCorrector method generally, so I have a little comment  idea  on it (and yes, I have no time until 21 Dec to check it): When you calculated the slope in the line 70 and you use the Euler estimation of Y, you get an estimated value of the slope in X+dX. But when Y new value is calculated, you have a possibility to improve this slope if you use the second order estimated Y value in X+dX. Maybe repeating the line 70 with the most precise Y value gives a better estimation?! Csaba 

12152019, 12:23 AM
Post: #11




RE: (71B) EulerTaylor method for the HP71B
Quote: Thanks. You're lucky man, you have more spare time than me. It is an endyearrush in power industry and today (Saturday) is a work day in Hungary (in this December two Saturdays are workdays  these are "relocated" to increase the Christmas holiday period). Repeating line 70 generates more error since that step has no theoretical justification. I did try it and the percent error jumped significantly!!! Namir 

12152019, 08:25 AM
(This post was last modified: 12152019 08:26 AM by Csaba Tizedes.)
Post: #12




RE: (71B) EulerTaylor method for the HP71B
(12152019 12:23 AM)Namir Wrote: Repeating line 70 generates more error since that step has no theoretical justification. I did try it and the percent error jumped significantly!!! Of course you must to store the temporary Euler estimation (YE) first: 50 D1=FND(X,Y) 60 X=X+H @ YE=Y+H*D1 Then use this, here the Y is totally same calculated as earlier: 70 Y=YE+H/2*(FND(X,YE)D1) We can repeat the calculation ("correction") of Y use the Euler estimation (YE) and the previously calculated, second order estimated Y: 75 Y=YE+H/2*(FND(X,Y)D1) Csaba 

12152019, 08:14 PM
Post: #13




RE: (71B) EulerTaylor method for the HP71B
(12152019 08:25 AM)Csaba Tizedes Wrote: 50 D1=FND(X,Y) The code looks good, except YE is an invalid variable. HP71B only allowed single character variable name, with optional single digits, like Y1 However, with 3 FND() calls per loop, it cost upto 50% more time than the original 2 FND() A more practical approach might be to keep original code (equivalent to removing line 75), but more points (smaller h) Another approach is to add "rough" 3rd order correction. Using backward differences, this only do 2 FND() calls per loop. h³ Y''' / 6 ≈ h³ (∇Y''/h) / 6 ≈ h² ∇(ΔY'/h) / 6 = ∇(½ h ΔY') / 3 Code: 50 D1=FND(X,Y) >RUN y = 1.64876111344 exact = 1.6487212707 %err = 2.4165843377E3 For comparison, using *exact* Y''' gives similar result Y' = X Y Y'' = X Y' + Y = X² Y + Y Y''' = (X² Y' + 2XY) + Y' = X³Y + 2XY + XY = XY(X² + 3) >75 D3=D1*((XH)^2+3) >77 Y=Y+H^3/6*D3 > >RUN y = 1.64876145117 exact = 1.6487212707 %err = 2.43706869767E3 3rd order correction patch shines when Y''' is relatively large. Example, with Y' = Y, X = 0 to 1 step 0.01 Exact Y = Y' = Y'' = Y''' = 2.71828182846 1st order correction, Y = 2.70481382938 2nd order correction, Y = 2.71823686266 3rd order correction, Y = 2.71828104622 

12162019, 11:23 PM
Post: #14




RE: (71B) EulerTaylor method for the HP71B
Thanks Albert for providing approximations for the seocnd and third derivative to extend the original Euler method.
If Albert ever attends an HHC conference then I am buying him dinner on Friday night (the eve of the start of the conference) OR Sunday night (where the HHC attendees go out for dinner). I am good for my word! Here is an update of my original listing to take into account Chan's suggestions. Code: 10 REM EULERTAYLOR WITH 3 DERIVATIVES 

« Next Oldest  Next Newest »

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