Calculators and numerical differentiation

[robve]

Came across this article that might be of interest to this forum: "calculators and numerical differentiation" http://blog.damnsoft.org/tag/fx-880p/

Some ways vintage calculators approximated differentiation "the right and wrong way".

- Rob

[Paul Dale]

The WP 34S doesn't use either of the methods discussed. It uses an order 10 method for it's derivative with fallbacks to order 6 and 4 if the function doesn't evaluate properly. For the second derivative is again uses an order 10 method with a fallback to an order 4 method. The order 10 methods are a weighted summation of \( f(x \pm 1) \), \( f(x \pm 2) \), \( f(x \pm 3) \), \( f(x \pm 4) \) and \( f(x \pm 5) \).

There are good reasons for not wanting to evaluate the function at the x specified.

[Albert Chan]

(10-30-2020 09:57 PM)robve Wrote:  Came across this article that might be of interest to this forum: "calculators and numerical differentiation" http://blog.damnsoft.org/tag/fx-880p/

The author suggested Casio is doing central difference formula, based on Casio CFX-9×50 manual.
On closer reading, it only *illustrated* what is central difference.

Using the example f(x)=1/x, a = 0.001, h = 0.0001

f'(a) ≈ (f(a+h) - f(a-h)) / (2h) = -1/(a²-h²) < -1/a²

Casio fx-570MS: d/dx(1/x, 0.001, 0.0001) = -999974.6848 > -1/a²
Casio fx-115ES+: d/dx(1/x, 0.001) = -999999.999994767 > -1/a²

This suggested Casio is not using central difference formula as-is.
Something more is involved ...

[Wes Loewer]

Instead of averaging the backwards and forwards differences, why not just check to see if f(a) exists and if it does then do the central difference method?

Notice in the TI manual that the default epsilon can be overridden. Same goes the their 8x numeric models.

The non-CAS Nspire apparently has a bit of CAS hidden under the hood because it does not use this approximation. It appears to evaluate the derivative symbolically and then evaluates that expression numerically, keeping the CAS carefully hidden from the user.

(Had anyone previously seen the del operator used for the backwards difference as shown in the Casio manual?)

[Albert Chan]

(11-01-2020 05:39 AM)Wes Loewer Wrote:  Instead of averaging the backwards and forwards differences, why not just check to see if f(a) exists and if it does then do the central difference method?

Try taken derivative of f(x) = |x|, at x=0.

f(0) = |0| = 0, thus exist.
Central difference slope = (|0+h| - |0-h|)/(2h) = 0

Forward difference slope = (h-0)/h = 1
Backward difference slope = (0-h)/h = -1
2 slopes does not match. f'(0) does not exist.

Oh ... Casio does not check this, d/dx(√(X²), 0) → 0

Quote:The non-CAS Nspire apparently has a bit of CAS hidden under the hood because it does not use this approximation. It appears to evaluate the derivative symbolically and then evaluates that expression numerically, keeping the CAS carefully hidden from the user.

How did you deduce there is hidden CAS under non-CAS Nspire ?
May be an example ?

Quote:Had anyone previously seen the del operator used for the backwards difference as shown in the Casio manual?

I think many text use the same symbols (or equivalent):
Ref: Fundamentals of Numerical Analysis, by Stephen Kellison

Code:

Operators:
Δf(x) = f(x+1) - f(x)       // forward difference
∇f(x) = f(x) - f(x-1)       // backward difference
Ef(x) = f(x+1)              // stepping
δf(x) = f(x+½) - f(x-½)     // central difference
µf(x) = (f(x+½)+f(x-½))/2   // average
Df(x) = f'(x)               // derivative

Or, operator notation: \(Δ=E-1,\quad ∇=1-E^{-1}\)

Combining µδ, we have the 3 points central difference slope formula

\( hD ≈ µδ = \large\left({E^{½} + E^{-½} \over 2}\right) \normalsize(E^{½}-E^{-½})
= \large{E\;-\;E^{-1} \over 2}
= \large{Δ\;+\;∇ \over 2}
\)

For more accuracy, we can add more terms: (note, there is no even powers of δ )

\( hD ≡ µδ \large\left(1 - {δ^2\over 6} + {δ^4\over 30} - {δ^6\over 140} + {δ^8\over 630}
- {δ^{10}\over 2772} + {δ^{12}\over 12012} - {δ^{14}\over 51480} \;+\; ... \right) \)

(10-30-2020 11:41 PM)Paul Dale Wrote:  The order 10 methods are a weighted summation of \( f(x \pm 1) \), \( f(x \pm 2) \), \( f(x \pm 3) \), \( f(x \pm 4) \) and \( f(x \pm 5) \).

Lets build D (order 10), without using central difference table.

XCas> c := E^(1/2) - E^(-1/2); // central difference
XCas> m := (E^(1/2) + E^(-1/2))/2; // mean
XCas> [mc, cc] := expand(simplify([m,c] .* c))       → [E/2 - 1/(E*2) , E - 2 + 1/E]

Note the symmetry of m*c, c*c. This suggested D also have similar symmetry

XCas> simplify(mc/h * horner([1/630, -1/140, 1/30, -1/6, 1], cc))

\(→ D = \large \frac {2100(E-E^{-1})\;
-\;600(E^2-E^{-2})\;
+\;150(E^3-E^{-3})\;
-\;25(E^4-E^{-4})\;
+\;2(E^5-E^{-5})}
{2520h}\)

[Albert Chan]

(11-01-2020 05:39 PM)Albert Chan Wrote:  For more accuracy, we can add more terms: (note, there is no even powers of δ )

\( hD ≡ µδ \large\left(1 - {δ^2\over 6} + {δ^4\over 30} - {δ^6\over 140} + {δ^8\over 630}
- {δ^{10}\over 2772} + {δ^{12}\over 12012} - {δ^{14}\over 51480} \;+\; ... \right) \)

We can show Df(0) = f'(0) does not require calculating f(0).
In other words, operator form will not have a constant term, the "1" operator.

From previous post, we have µδ = (E-1/E)/2, δδ = (E+1/E) - 2
Doing "operator" mathematics, with x = log(E), we have:

µδ = sinh(x)
δδ = 2*cosh(x) - 2

Hyperbolics identities:
(1): cosh(z1)*cosh(z2) = (cosh(z1 - z2) + cosh(z1 + z2)) / 2
(2): sinh(z1)*cosh(z2) = (sinh(z1 - z2) + sinh(z1 + z2)) / 2

hD = sinh(x) * (k1 + k2*cosh(x) + k3*cosh(x)^2 + k4*cosh(x)^3 + ...)       // apply (1)
     = sinh(x) * (k1' + k2'*cosh(x) + k3'*cosh(2x) + k4'*cosh(3x) + ... )       // apply (2)
     = k1''*sinh(x) + k2''*sinh(2x) + k3''*sinh(3x) + k4''*sinh(4x) + ...

sinh(nx) = (Eⁿ - E^-n)/2

→ this explained why Eⁿ coefs = negative of E^-n coefs.
→ RHS terms will not generate constant term (i.e., no "1" operator)

→ D does not require calculating f(0)
→ same for D^odd_powers, since RHS is still linear combinations of sinh's.

[Wes Loewer]

(11-01-2020 05:39 PM)Albert Chan Wrote:  How did you deduce there is hidden CAS under non-CAS Nspire ?
May be an example ?

I may have made a leap in my logic, but the fact that examples such as |x| or 1/x at x=0 or at x=0.0001 produce the correct answer in the non-CAS Nspire while these are incorrect on the 84+ lead me to believe that the non-CAS model must be doing something CAS-like under the hood. I had never come across a counter-example. It also made sense to me the it would be very easy to share the same code as the Nspire-CAS for such calculations.

You prompted me to look in the Nspire manual which gives some insight.

Quote:nDerivative()
Returns the numerical derivative calculated using auto differentiation methods.
...
Note: The nDerivative() algorithm has a limitiation [sic]: it works recursively through the unsimplified expression, computing the numeric value of the first derivative (and second, if applicable) and the evaluation of each subexpression, which may lead to an unexpected result.

Consider the example on the right. The first derivative of x•(x^2+x)^(1/3) at x=0 is equal to 0. However, because the first derivative of the subexpression (x^2+x)^(1/3) is undefined at x=0, and this value is used to calculate the derivative of the total expression, nDerivative() reports the result as undefined and displays a warning message.

If you encounter this limitation, verify the solution graphically. You can also try using centralDiff().

(The Npsire has a centralDiff() function as well that behaves like the 84+.)

The part starting with "Note: " is not in the Nspire CAS manual. The Nspire CAS gives the correct answer for this example.

When I first saw the above example, I thought that I must have been wrong about the numeric model having an internal CAS since the numeric model does not give the correct answer while the CAS model does. However, now that I'm reading it again, I'm thinking that I may have been right after all. The fact that it says "the subexpression (x^2+x)^(1/3) is undefined at x=0" means that the calculator must be breaking the expression down into subexpressions and evaluating their derivatives (consistent with the product rule) which means that it must have some CAS capabilities rather than just evaluating the whole expression numerically.

So my current thinking is that the Nspire must have at least some CAS capabilities under the hood, but not to the extent as the Nspire CAS.

Thoughts?

[CMarangon]

Hello!

Nice post!

I suggest you to increase the font size, because we have many who access
our pages with cell phones or other guys like me who don't like to wear glasses.:-)

[Albert Chan]

Slightly off topics, for f(x) = x*g(x), getting f'(0) is easier taking limit directly.

\(f(x) = x·g(x) = x·\sqrt[3]{x^2+x}\)

\(f'(0) = \displaystyle{\lim_{h \to 0}} {f(h)-f(0)\over h}
= \displaystyle{\lim_{h \to 0}}\; g(h) = g(0) = 0\)

Getting derivative take more work, and easier to make mistakes.

f' = (x*g)' = g + x*g'

g has the form z^n, where z = x²+x, n = 1/3

g' = (n*z^(n-1)) * z' = (n*g) * (z'/z) = g/3 * (2x+1)/(x²+x)

f' = g + g/3 * (2x+1)/(x+1) = g * (5x+4)/(3x+3)

f'(0) = g(0) * 4/3 = 0

---

Another way, by shape of the curve.

\(f(x) = x·\sqrt[3]{x^2+x} = \large\frac{\sqrt[3]{(x^2+x)^4}}{x+1}\)

f(ε) > 0, f(-ε) > 0, f(0) = 0

→ f(0) is local minimum
→ f'(0) = 0

[Wes Loewer]

(11-03-2020 06:55 PM)CMarangon Wrote:  I suggest you to increase the font size, because we have many who access our pages with cell phones or other guys like me who don't like to wear glasses.:-)

I'm not sure what you are seeing. My post and all the other posts have the same font size. I checked on my phone as well and they all look the same.

[Wes Loewer]

(11-03-2020 10:14 PM)Albert Chan Wrote:  Slightly off topics, for f(x) = x*g(x), getting f'(0) is easier taking limit directly.

\(f(x) = x·g(x) = x·\sqrt[3]{x^2+x}\)

\(f'(0) = \displaystyle{\lim_{h \to 0}} {f(h)-f(0)\over h}
= \displaystyle{\lim_{h \to 0}}\; g(h) = g(0) = 0\)

If you pull out an \(x\), then \( x \cdot (x^2+x)^{1/3}\) becomes \( x^{4/3} \cdot (x+1)^{1/3}\) which the non-CAS Npsire handles correctly.