(71B) Euler-Taylor method for the HP-71B

[Namir]

Euler-Taylor 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)

HP-71B Listing
==============

Code:

10 REM EULER METHOD WITH 2 DERIVATIVES
20 DEF FNE(X)=EXP(0.5*x^2)
30 DEF FNX(X,Y) = X*y
40 DEF FND(X,Y)=X*FNX(X,Y)+y
50 READ A,B,Y,H
60 DATA 0,1,1,0.01
70 N=INT((B-A)/H+.5)
80 X=A
90 DISP "WAIT PLEASE..." @ WAIT .1
100 FOR I=1 TO N
110 Y = Y + H*FNX(X,Y) + H^2/2*FND(X,Y)
120 X = X + H
130 NEXT I
140 DISP "Y=";Y @ PAUSE
150 Y0 = FNE(B)
160 DISP "Y EXACT=";Y0 @ PAUSE
170 DISP "%ERR=";100*(Y-Y0)/Y0
180 END

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


HP-71B 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
20 DEF FNE(X)=EXP(0.5*x^2)
30 READ A,B,Y,H
40 DATA 0,1,1,0.01
50 N=INT((B-A)/H+.5)
60 X=A
70 DISP "WAIT PLEASE..." @ WAIT .1
80 FOR I=1 TO N
90 D1 = X*Y
100 D2 = X*D1+Y
110 Y = Y + H*D1 + H^2/2*D2
120 X = X + H
130 NEXT I
140 DISP "Y=";Y @ PAUSE
150 Y0 = FNE(B)
160 DISP "Y EXACT=";Y0 @ PAUSE
170 DISP "%ERR=";100*(Y-Y0)/Y0
180 END

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 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

[Csaba Tizedes]

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

[Namir]

(12-10-2019 06:02 PM)Csaba Tizedes Wrote:  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

1) Thatls for spotting the error. I corrected it in both listings.
2) HP-71B 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

[Albert Chan]

(12-11-2019 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((X1-X)/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*(Y-Y1)/Y1

>RUN
y       = 1.64871423741
exact = 1.6487212707
%err =-4.26590602365E-4

Edit:To compare apples to apples, I removed the assumption final X=X1
So, both exact and estimate uses final X as end-point.

line 90: Y1=EXP(.5*X1*X1) changed to Y1=EXP(.5*X*X)

[Namir]

Thank you very much Albert for your approximation and the updated HP-71B listing. My hats off to you sir!!

:-)

Namir

[Namir]

(12-13-2019 02:12 PM)Albert Chan Wrote:  

(12-11-2019 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 ≈ 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((X1-X)/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*X1*X1)
100 DISP "y =";Y
110 DISP "exact=";Y1
120 DISP "%err =";100*(Y-Y1)/Y1

>RUN
y       = 1.64871423741
exact = 1.6487212707
%err =-4.26590602365E-4

Do I have your permission to update/adapt my HP-71B listing using your approximation for the second derivative? I will give you credit for that update (and also for the update of the HP-41C Listing).

Namir

[Albert Chan]

(12-13-2019 04:10 PM)Namir Wrote:  Do I have your permission to update/adapt my HP-71B 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.266e-4 to to 0.6612

2nd order correction improved accuracy 1550X

[Namir]

Albert,

For last last month I have been studying numerical methods for ODE that are based on Runge-Kutta 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 Runge-Kutta 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

[Namir]

Adapting my original HP-71B to Albert Chan's approximaton of the second derivative, we get:

Code:

10 REM EULER METHOD WITH 2 DERIVATIVES
20 DEF FNE(X)=EXP(0.5*x^2)
30 DEF FNX(X,Y) = X*Y
40 READ A,B,Y,H
50 DATA 0,1,1,0.01
60 N=INT((B-A)/H+.5)
70 X=A
80 DISP "WAIT PLEASE..." @ WAIT .1
90 FOR I=1 TO N
100 D1=FNX(X,Y)
110 X = X+H @ Y= Y + H*D1
120 Y = Y + H/2*(FNX(X,Y)-D1)
130 NEXT I
140 DISP "Y=";Y @ PAUSE
150 Y0 = FNE(B)
160 DISP "Y EXACT=";Y0 @ PAUSE
170 DISP "%ERR=";100*(Y-Y0)/Y0
180 END

[Csaba Tizedes]

(12-13-2019 02:12 PM)Albert Chan Wrote:  Δy' ≈ y'(x+h, y + h y'(x,y)) - y'(x,y)
→ y ≈ y + h y' + ½ h² y'' ≈ y + h y' + ½ h Δy'

10 DEF FND(X,Y)=X*Y
...
50 D1=FND(X,Y)
60 X=X+H @ Y=Y+H*D1
70 Y=Y+H/2*(FND(X,Y)-D1)
...

Thanks. You're lucky man, you have more spare time than me. It is an end-year-rush 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 Predictor-Corrector 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

[Namir]

Quote: Thanks. You're lucky man, you have more spare time than me. It is an end-year-rush 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 Predictor-Corrector 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

Repeating line 70 generates more error since that step has no theoretical justification. I did try it and the percent error jumped significantly!!!

Namir

[Csaba Tizedes]

(12-15-2019 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

[Albert Chan]

(12-15-2019 08:25 AM)Csaba Tizedes Wrote:  50 D1=FND(X,Y)
60 X=X+H @ YE=Y+H*D1
70 Y=YE+H/2*(FND(X,YE)-D1)
75 Y=YE+H/2*(FND(X,Y)-D1)

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)
60 X=X+H @ Y=Y+H*D1
70 D2=H/2*(FND(X,Y)-D1) @ Y=Y+D2
75 IF I>1 THEN Y=Y+(D2-D3)/3
77 D3=D2

>RUN
y       = 1.64876111344
exact = 1.6487212707
%err = 2.4165843377E-3

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*((X-H)^2+3)
>77 Y=Y+H^3/6*D3
>
>RUN
y       = 1.64876145117
exact = 1.6487212707
%err = 2.43706869767E-3

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

[Namir]

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 EULER-TAYLOR WITH 3 DERIVATIVES
20 DEF FNX(X,Y)=X*Y
30 FNE(X) = EXP(0.5*X^2)
40 READ A, B, Y, H
50 DATA 0, 1, 1, 0.01
60 N=INT((X1-X)/H+.5)
70 X=A
80 FOR I=1 TO N
90 D1=FNX(X,Y)
100 X=X+H @ Y=Y+H*D1
110 D2=H/2*(FNX(X,Y)-D1) @ Y=Y+D2
120 IF I>1 THEN Y=Y+(D2-D3)/3
130 D3=D2
140 NEXT I
150 Y1=FNE(B)
160 DISP "Y =";Y
170 DISP "Y EXACT=";Y1
180 DISP "%ERR =";100*(Y-Y1)/Y1
190 END