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

[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