(15C) (DM 15L) Analyzing a function f(x)

[rawi]

This program helps in analyzing a function f(x) in a given area of x. The function should be defined and differentiable in the whole area we want to analyze.
It does the following tasks:
- Make a table of function values for sketching the curve
- Searching for zero points
- Searching for extrema
- Searching for inflection points (extrema of the first derivative)
- Computing the function value, the 1st and 2nd derivative at any given value for x
Method:
The program approximates the derivate by adding and subtracting the small value d stored in Register I to x and evaluation the function. The derivative is estimated by (f(x+d)-f(x-d)/(2d).
The second derivative is approximated in a similar way by using additionally f(x).
Making the small value too big or too small delivers poor approximations. In this example 0.01 delivers good results with four significant digits after decimal point.
Use of registers: Register 0-3 and I are used.

How to apply:

Preparation:
- Put in the function under LBL E. It has to take the x-value from the X-register and return f(x) in X register.
- Put in a small number (e.g. 0.01) in Register I. This is used for the numerical computation of the first and second derivates.

Application:
1. To generate a table of function values:
- Starting point of x ENTER
- End point of X ENTER
- Step length x-values
- Press f A
Output:
- f(x) at starting point R/S
- f(x) at starting point + 1 * step length R/S
- f(x) at starting point + 2* step length R/S
….
Note: The Y-register has X, the X-register f(x)
So at every step you can check the value of x by pressing x<>y

2. To find a zero point of the function: Put in two estimates in the stack and press f SOLVE E (OK, this was easy)
3. To find an extremum: Put in two estimates in stack and press f SOLVE B. Program returns x value of extremum. To get f(x) at extremum press f E.
4. To find an inflection point (i.e. an extremum of the first derivate) put in two estimates in stack and press f SOLVE C. Program returns x value of inflection point. To get f(x) at this point press f E. If no inflection point is found program returns “Error 8”.
5. To get the first derivatie at a given x: Put in x in stack, press f B.
6. To get the second derivative of a given x: Put in x in stack, press f C.

Example:
Let’s look at the function f(x) = (e^(-x²)+x²-5)/(sin(x)+2) in the area from -4 to +4.

Preparation:
The program for LBL E is listed below. We put 0.01 in STO I.

1) We want to generate a table from x=-4 to x=4 with step width 1:
-4 ENTER 4 ENTER 1 f A: 3.9901 (This is the value of f(x) at -4. To check x press x<>y: -4)
R/S 2.1519 (=f(-3))
R/S -0.9000
R/S -3.1351
R/S -2.0000
R/S -1.2783
R/S -0.3374
R/S 1.8682
R/S 8.8482, which equals f(4).

2. From the table we know that there has to be a zero point between x=-3 and x=-2. Let’s search the zero-point: Put in stack: -3 ENTER -2 f SOLVE E. Result: -2.2346 (x-value of zero point).
There is another zero-point between x=2 and x=3: 2 ENTER 3 f SOLVE E Result: 2,2346

3. From the table we know that there is a minimum between x=-2 and x=0. Let’s search for the extremum:
-2 ENTER 0 f SOLVE B Result: -1.0795. To get f(x) at minimum press f E -> -3.1503

4. There can be an inflection points between -4 and -3, -3 and -2, -1 and 0 and 1 and 3. Let’s check that:
-4 ENTER -3 f SOLVE C -> -3,5395. For f(x) at inflection point press f E -> 3.1531
-3 ENTER -2 f SOLVE C -> -1.9835. f E -> -0.9680
-1 Enter 0 f SOLVE C -> -0.4898 f E -> -2.5978
0 ENTER 3 f SOLVE C -> 0.6831 f E -> -1.4846

Enjoy!
Raimund


01: f LBL A
02: STO 0
03: R↓
04: STO 2
05: x<>y
06: STO 1
07: f LBL 0
08: GSB E
09: RCL 1
10: x<>y
11: R/S
12: RCL 2
13: RCL 1
14: RCL 0
15: +
16: STO 1
17: x<=y?
18: GTO 0
19: g RTN
20: f LBL B
21: GSB 1
22: -
23: RCL I
24: 2
25: X
26: ./.
27: g RTN
28: f LBL C
29: STO 1
30: GSB E
31: STO 0
32: RCL 1
33: GSB 1
34: RCL 0
35: -
36: RCL I
37: ./.
38: x<>y
39: RCL 0
40: x<>y
41: -
42: RCL I
43: ./.
44: -
45: RCL I
46: ./.
47: g RTN
48: f LBL 1
49: STO 2
50: RCL I
51: +
52: GSB E
53: STO 3
54: RCL 2
55: RCL I
56: -
57: GSB E
58: RCL 3
59: x<>y
60: g RTN
61: f LBL E
62: ENTER
63: ENTER
64: ENTER
65: g x²
66: CHS
67: e^x
68: x<>y
69: g x²
70: +
71: 5
72: -
73: x<>y
74: SIN
75: 2
76: +
77: ./.
78: g RTN

[Valentin Albillo]

.
Hi, rawi:

First of all, welcome to the MoHPC fora. I'm sure you'll enjoy your stay here.

(06-05-2020 05:04 PM)rawi Wrote:  This program helps in analyzing a function f(x) in a given area of x.

Thanks for posting your efforts, that's the spirit. Some comments:

Quote:The program approximates the derivate by adding and subtracting the small value d stored in Register I to x and evaluation the function. The derivative is estimated by (f(x+d)-f(x-d)/(2d).
The second derivative is approximated in a similar way by using additionally f(x).
Making the small value too big or too small delivers poor approximations. In this example 0.01 delivers good results with four significant digits after decimal point.

First, making the increment equal to 0.01, a constant, for all x arguments isn't the best idea, it must really depend on the value of x, i.e.: if you're computing the derivative at x=5 then the increment might be 5*0.01 = 0.05 so you compute the derivative using f(5+0.05) and f(5-0.05). If the argument is x=100, you'd compute f(100+1) and f(100-1), and so on. For x=0, simply use 0.01.

Second, the HP-15C is a 10-digit machine so you'll get greater accuracy using a constant smaller than 0.01. For the first derivative, use 0.00001 instead, and for the second derivative use 0.001. The first is 10^Int(-10[digits]/2) while the second is 10^Int(-10[digits]/3).

Third, for the HP-15C specifically there are better ways to compute first derivatives, faster and giving up to 10 correct digits. See how it's done by having a look at this PDF document:

Technique for computing first order derivatives using complex operations

Best regards and have a nice weekend.
V.

[rawi]

Hi Valentino,
first of all thank you very much for your warm welcome and your very helpful remarks.

I tried to implement the result of the challenge but it did not give the right results. Perhaps this is because I have a DM15L and not the original HP15C?

Concerning the amount of the small value d: You will find a file attached with the results for the first and the second derivative for two values of my function and as well for the three examples of the challenge with d ranging from 1E-6 to 1E-2. The result: In all cases the result was best for the second derivative for d=0.01. For the first derivative it was better when it was smaller than that (as you said). It was best for d=0.001 or d=0.0001. I also checked d=0.1 which had always worse results as well for the first as the second derivative.

So I changed my code that for the derivative the d is divided by 10.

You mentioned that the difference should change if x changes. I am not sure, whether this is correct for all functions. Think of f(x)=(sin(x))^3+(cos(x))^3, which is always in the same range regardless whether x is zero or 100. Why should d be changed?

Thank you very much once more and

Best regards

Raimund

[Valentin Albillo]

.
Hi again:

(06-06-2020 11:10 PM)rawi Wrote:  Hi Valentino,

Bad start. It's "Valentin", not "Valentino". I'm not Italian, you know.

Quote:first of all thank you very much for your warm welcome and your very helpful remarks.

You're welcome.

Quote:I tried to implement the result of the challenge but it did not give the right results. Perhaps this is because I have a DM15L and not the original HP15C?

As far as I know, the SwissMicros' DM15 and DM15L use a modified original HP-15C ROM so they must work Ok and produce the same results as in my challenge. Perhaps you're doing something wrong (RAD mode ? Complex mode ?)

Quote:Concerning the amount of the small value d: You will find a file attached with the results [...]

Sorry but you attached a .docx (i.e., MS Word) document and I'm reading the forum using an Android tablet so I can't open such documents, attach a PDF document instead.

Quote:The result: In all cases the result was best for the second derivative for d=0.01. For the first derivative it was better when it was smaller than that (as you said). It was best for d=0.001 or d=0.0001. I also checked d=0.1 which had always worse results as well for the first as the second derivative. So I changed my code that for the derivative the d is divided by 10.

Good. I conducted some tests, like these:

f(x)=sin(x), 12-digit machine:

      inc       x       f'(x)       cos(x)       err
      ----------------------------------------------------------------------
      0.1       1 0.5394022522 0.5403023059 -0.0009000537
      0.01      1 0.5402933009 0.5403023059 -0.0000090050
      0.001     1 0.5403022158 0.5403023059 -0.0000000901
      0.0001    1 0.5403023050 0.5403023059 -0.0000000009
      0.00001   1 0.5403023059 0.5403023059  0.0000000000 <<<

f(x)=sin(x), 10-digit machine:

      inc       x       f'(x)       cos(x)       err
      ----------------------------------------------------------------------
      0.1       1 0.539402253 0.540302306 -0.000900053
      0.01      1 0.540293300 0.540302306 -0.000009006
      0.001     1 0.540302200 0.540302306 -0.000000106  <<<
      0.0001    1 0.540302000 0.540302306 -0.000000306
      0.00001   1 0.540305000 0.540302306  0.000002694

f"(x) = (f(x+h)-2*f(x)+f(x-h))/h^2 f(x)=sin(x), f"(x)=-sin(x), 10-digit machine

      inc       x       f"(x)       -sin(x)       err
      ----------------------------------------------------------------------
      0.1       1 -0.840770000 -0.841470985 0.000700985
      0.01      1 -0.841470000 -0.841470985 0.000000985 <<<
      0.001     1 -0.842000000 -0.841470985 0.000529015

f(x)=e^x, f"(x)=e^x, 10-digit machine

      0.1       1 1.106092200 2.718281828 -1.612189628
      0.01      1 2.718310000 2.718281828  0.000028172 <<<
      0.001     1 2.719000000 2.718281828  0.000718172

which essentially agree with your observations.

On another note, your program can be somewhat optimized in various places as it seems you're not using some of the advanced capabilities of the HP-15C, as for example recall arithmetic. For instance, your code computing the second derivative (LBL C, steps 28-47, 20 steps in all) can be shortened to just 12 steps (i.e. 40% shorter) like this:

28: f LBL C   stack contents (X Y)
29: STO 1     x
30: GSB E     f(x)
31: STO 0     f(x)
32: RCL 1     x
33: GSB 1     f(x-inc)     f(x+inc)
34: +         f(x-inc) + f(x+inc)
35: RCL- 0    f(x-inc) + f(x+inc) - f(x)
36: RCL- 0    f(x-inc) + f(x+inc) - 2*f(x)
37: RCL/ I    [f(x-inc) + f(x+inc) - 2*f(x)]/inc
38: RCL/ I    [f(x-inc) + f(x+inc) - 2*f(x)]/inc^2
39: g RTN

Have a nice week.

V.

[rawi]

Hi Valentin,

First my apologies for misspelling your name.

I replaced the word-cocument by a pdf document. And I put in a new version of the program
which is much improved due to your valuable tips in a second post. Thank you very much for that.

Best

Raimund

[rawi]

This is an improved version of the program posted at the beginning of this thread.

This program helps in analyzing a function f(x) in a given area of x. The function should be defined and differentiable in the whole area. The program does the following tasks:
- Make a table of function values for sketching the curve
- Search for zero points
- Search for extrema
- Search for inflection points
- Compute the function value, the 1st and 2nd derivative at any given value for x
- Compute the length of the function curve between to given x-values

The improvements are:
- For the computation of the first and the second derivative different increment values are used resulting in more precise results.
- Computation of function length is included.
- Code was shortened so that the length of the program did not increase much despite the improvements.
Many of these improvements are due to suggestions by Valentin Albillo whom I want to thank here.

Method:
The program approximates the derivative by adding and subtracting the small value d stored in Register I to x and evaluation the function. The derivative is estimated by (f(x+d)-f(x-d)/(2d).
The second derivative is approximated in a similar way by using additionally f(x) and a 10 times greater d to get numerical stable results.
Making the small value too big or too small delivers poor approximations. It is recommended to set it by .001. By experience the deviation then is about 1 unit in the fifth figure.
Use of registers: Registers 0-3 and I are used.

How to apply:
Preparation:
- Put in the function under LBL E. It must take the x-value from the X-register and return f(x) in X register.
- Put in a small number (recommendation: 0.001) in Register I.
Application:
1. To generate a table of function values:
- Starting point of x ENTER
- Ending point of X ENTER
- Step length x-values
- Press f A
Output:
- f(x) at starting point R/S
- f(x) at starting point + 1 * step length R/S
- f(x) at starting point + 2* step length R/S
….
Note: In the Y-register there is the x-value, in the X-register the value of f(x)
So at every step you can check the value of x by pressing x<>y
2. To find a zero point of the function: Put in two estimates in the stack and press f SOLVE E (OK, this was easy)
3. To find an extremum: Put in two estimates in stack and press f SOLVE B. Program returns x value of extremum. To get f(x) at extremum press f E. If no extremum is found program returns “Error 8”.
4. To find an inflection point (i.e. an extremum of the first derivate) put in two estimates in stack and press f SOLVE C. Program returns x value of inflection point. To get f(x) at this point press f E. If no inflection point is found program returns “Error 8”.
5. To get the first derivative at a given x: Put in x in stack, press f B.
6. To get the second derivative of a given x: Put in x in stack, press f C.
7. To get the curve length of the function between two points: Put lower x in Y-register and upper x in X register and press integration key D.

Example:
Function f(x) = (e^(-x²)+x²-5)/(sin(x)+2) in the area from -4 to 4.

Preparation:
The program for LBL E is listed below. 0.001 STO I.
1) Generate a table from x=-4 to x=4 with step width 1:
-4 ENTER 4 ENTER 1
f A: 3.9901 (This is the value of f(x) at -4. To check x press x<>y: -4)
R/S 2.1519 (=f(-3))
R/S -0.9000
R/S -3.1351
R/S -2.0000
R/S -1.2783
R/S -0.3374
R/S 1.8682
R/S 8.8482, which equals f(4).
2) Search for zero point between x=-3 and x=-2. -3 ENTER -2 f SOLVE E Result: -2.2346 (x-value of zero point). Search for zero point between x=2 and x=3: 2 ENTER 3 f SOLVE E -> 2,2346
3) Search for minimum between x=-2 and x=0. -2 ENTER 0 f SOLVE B -> -1.0796. To get f(x) at minimum press f E -> -3.1503
4) Search for inflection points:
-4 ENTER -3 f SOLVE C -> -3,5395. For f(x) at inflection point press f E -> 3.1531
-3 ENTER -1 f SOLVE C -> -1.9825. f E -> -0.9690
-1 ENTER 0 f SOLVE C -> -0.4898 f E -> -2.5978
0 ENTER 3 f SOLVE C -> 0.6831 f E -> -1.4846
5) Estimation the length of the function curve from -4 to 4:
-4 ENTER 4 f INTEGRATION KEY D -> 21.5315

PHP Code:

f LBL A        001 – 42,21,11
STO 0        002 – 44 0
R↓        003 – 33 
STO 2        004 – 44 2
x<>y        005 – 34 
STO 1        006 – 44 1
f LBL 0        007 – 42,21,0 
GSB E        008 – 32 15
RCL 1        009 – 45 1
x<>y        010 – 34 
R/S        011 – 31 
RCL 0        012 – 45 0
STO+1        013 – 44,40,1
RCL 2        014 – 45 2
RCL 1        015 – 45 1
g x<=y        016 – 43 10
GTO 0        017 – 22 0
g RTN        018 – 43 32
f LBL B        019 – 42,21,12         
GSB 1        020 – 32 1
-        021 – 30 
RCL I        022 – 45 25
/        023 – 10 
2        024 – 2 
/        025 – 10  
g RTN        026 – 43 32
f LBL C        027 – 42,21,13
EEX        028 – 26 
1        029 – 1 
STOXI        030 – 44,20,25
R↓         031 – 33 
STO 1        032 – 44 1
GSB E        033 – 32 15
CHS        034 – 16 
STO 0        035 – 44 0
STO+0        036 – 44,40,0
RCL 1         037 – 45 1
GSB 1        038 – 32 1
+        039 – 40
STO+0        040 – 44,40,0 
RCL I        041 – 45 25
STO/0        042 – 44 10,0 
STO/0        043 – 44,10,0
EEX        044 – 26 
1        045 – 1
STO/I        046 – 44,10,25 
RCL 0        047 – 45 0
g RTN        048 – 43 32
f LBL 1        049 – 42,21,1
STO 2        050 – 44 2
RCL I        051 – 45 25
-        052 – 30 
GSB E        053 – 32 15
STO 3        054 – 44 3
RCL 2        055 – 45 2
RCL I        056 – 45 25
+        057 – 40 
GSB E        058 – 32 15
RCL 3        059 – 45 3
g RTN        060 – 43 32
f LBL D        061 – 42,21,14
GSB B        062 – 32 12
g x²        063 – 43 11
1        064 – 1 
+        065 – 40 
SQRT        066 – 11 
g RTN        067 – 43 32
f LBL E        068 – 42,21,15
ENTER        069 – 36 
ENTER        070 – 36 
ENTER        071 – 36 
g x²        072 – 43 11
CHS        073 – 16
e^x        074 – 12 
x<>y        075 – 34 
g x²        076 – 43 11
+        077 – 40 
5        078 – 5 
-        079 – 30 
x<>y        080 – 34 
SIN        081 – 23 
2        082 – 2 
+        083 – 40 
/        084 – 10 
g RTN        085 – 43 32