Second derivative with complex numbers - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html) +--- Forum: General Forum (/forum-4.html) +--- Thread: Second derivative with complex numbers (/thread-16299.html) Second derivative with complex numbers - peacecalc - 02-09-2021 06:06 AM Hello all, I'm shure, in former threads has discussed that topic. After programming the first derivative with complex numbers in my DM42, I tried the same with the second derivative. For rememberance, the taylor-series for the first derivative (truncating after the first derivative): $f(x+i\cdot h) \approx f(x) + f'(\tau)|_{\tau = x}\cdot i h \quad \hbox{ take the imaginary part with}\quad \Im(f(x)) = 0 \quad \hbox{ because it is only real, one gets:} \quad f'(x) \approx \frac{1}{h}\Im(f(x+ih))$ Now let us see the taylor series until the second derivative (and then truncating further terms): $f(x+i\cdot h) \approx f(x) + f'(\tau)|_{\tau = x}\cdot i h + \frac{1}{2} f''(\tau)|_{\tau = x}\cdot (i h)^2 \quad \hbox{with} \quad i^2 = - 1 \Rightarrow f(x+i\cdot h) \approx f(x) + f'(\tau)|_{\tau = x}\cdot i h - \frac{1}{2} f''(\tau)|_{\tau = x}\cdot h^2$ now we take the real part: $\Re(f(x+i\cdot h)) \approx f(x) - \frac{1}{2} f''(\tau)|_{\tau = x}\cdot h^2 \quad\hbox{the term:} \quad f'(\tau)|_{\tau = x}\cdot i h \quad \hbox{vanish because it is only imagenary and one gets:} \quad f''(x) \approx - \frac{2}{h^2} \cdot \left( \Re(f(x+i\cdot h)) - f(x) \right)$ Maybe it is an advantage, because you only have to calculate one difference instead of three if you calculate the first order derivative without doing that in the complex plane. Please pay attention to, only the real part is taken from f(x+ih). Sincerely peacecalc RE: Second derivative with complex numbers - Albert Chan - 02-09-2021 03:45 PM Here is another way to obtain the derivative formulas $$\displaystyle f'(x) ≈ {f(x+ih) - f(x) \over ih}$$ $$\displaystyle f''(x) ≈ {f'(x) - f'(x-ih) \over ih} ≈ \frac{{{f(x+ih) - f(x) \over ih}} - {{f(x) - f(x-ih) \over ih}}}{ih} = {f(x+ih) - 2 f(x) + f(x-ih) \over -h^2}$$ Weierstrass approximation theorem: we can assume f(x) is polynomial. Conjugate of polynomial is polynomial of conjugate: f(x-hi) = conj(f(x+hi)) If f(x) is smooth and real, so does its derivatives → RHS must be real. $$\displaystyle f'(x) ≈ {\Im(f(x+ih)) \over h}$$ $$\displaystyle f''(x) ≈ {\Re(f(x+ih) - f(x)) \over -h^2/2}$$ --- Using central difference derivative formula, we also get the same f'(x) formula. Note: below does not assume real f'(x), RHS imaginary parts really cancelled out. $$\displaystyle f'(x) ≈ {f(x+ih) - f(x-ih) \over 2ih} = {\Im(f(x+ih)) - \Im(f(x-ih)) \over 2h} = {\Im(f(x+ih)) \over h}$$ → both derivative formulas should give accuracy similar to central difference formulas. XCas> f(x) := x*(x^3+5*x^2-21*x)﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ // example from here XCas> f'(4)﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 328 XCas> (f(x+h) - f(x-h)) / (2h) | x=4, h=1e-3 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 328.000021 XCas> im(f(x+h*i)) / h ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ | x=4, h=1e-3 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿→ 327.999979 XCas> f''(4)﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 270 XCas> (f(x+h)-2*f(x)+f(x-h))/h^2 | x=4, h=1e-3﻿ ﻿ → 270.000002132 XCas> re(f(x+h*i)-f(x)) / (-h^2/2) | x=4, h=1e-3 ﻿ → 269.999998125 RE: Second derivative with complex numbers - peacecalc - 02-11-2021 09:49 AM Hello Albert, thank you for your lucide comments about that topic. In deed I always thought that the way into the complex plane is an advantage in precision, but it seems not. What a pity! RE: Second derivative with complex numbers - Werner - 02-11-2021 12:36 PM It may avoid roundoff if you want high accuracy. Try the above examples on a real 42S, with h=1e-11. The complex first derivative gets 328 exactly, the difference function gets 300. Cheers, Werner RE: Second derivative with complex numbers - Albert Chan - 02-11-2021 02:13 PM (02-11-2021 12:36 PM)Werner Wrote:  It may avoid roundoff if you want high accuracy. Yes, but only with first derivative. (assumed im(f(x+hi)) does not hit with catastrophic cancellation) For second derivative, you still hit with subtraction cancellation errors. XCas> f(x) := x*(x^3+5*x^2-21*x) XCas> f1(x,h) := re((f(x+h)-f(x-h))/(2h)) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ // central difference 1st derivative XCas> f'(x) .- [f1(x,h), f1(x,h*i)] | x=4, h=1e-3 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → [-2.10000508787e-05, 2.10000000607e-05] XCas> f'(x) .- [f1(x,h), f1(x,h*i)] | x=4, h=1e-6 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → [-1.0320036381e-08, 2.09752215596e-11] XCas> f2(x,h) := re((f(x+h)-2*f(x)+f(x-h))/h^2) ﻿ ﻿ ﻿ ﻿ // central difference 2nd derivative XCas> f''(x) .- [f2(x,h), f2(x,h*i)] | x=4, h=1e-3 ﻿ ﻿ ﻿ ﻿ ﻿ → [-2.13206476474e-06, 1.87539626495e-06] XCas> f''(x) .- [f2(x,h), f2(x,h*i)] | x=4, h=1e-6 ﻿ ﻿ ﻿ ﻿ ﻿ → [0.0221821205923, 0.0506038300227] --- If |h| is not too small, we avoided catastrophic cancellation, error = O(h^2) We may take advantage of it, with h = ε*√i, which make h^2 purely imaginary. XCas> f'(x) - f1(x,h) | x=4, h=1e-3*√i ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → -2.44426701101e-12 XCas> f''(x) - f2(x,h) | x=4, h=1e-3*√i ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → -4.78621586808e-11 RE: Second derivative with complex numbers - Werner - 02-11-2021 05:49 PM You never cease to amaze me, Albert! One more for my files, complex derivatives with h^2 purely imaginary.. Werner RE: Second derivative with complex numbers - Werner - 02-11-2021 07:01 PM A Free42 program for both first and second derivative. I took h=sqrt(i), that seems to work as well, and then the division by h^2 for the second derivative is no longer necessary ;-) Update: whenever I think I'm being smart, it turns out I make things worse. So the program below uses h = sqrt(i)/1e3.. Update 2 : shorter Instructions: put the name of the function to derive in the ALPHA reg, X-value in X, press either F' or F" for first and second derivative, respectively. Code: 00 { 86-Byte Prgm } 01▸LBL 02 02 RCL+ "X" 03 ASTO ST L 04 GTO IND ST L 05▸LBL "F'" 06 AOFF 07 GTO 00 08▸LBL "F"" 09 AON 10▸LBL 00 11 LSTO "X" 12 -1ᴇ-12 13 SQRT 14 SQRT 15 LSTO "h" 16 XEQ 02 17 LSTO "." 18 RCL "h" 19 +/- 20 XEQ 02 21 FC? 48 22 GTO 00 23 STO+ "." 24 CLX 25 XEQ 02 26 STO+ ST X 27▸LBL 00 28 +/- 29 RCL+ "." 30 RCL "h" 31 FC? 48 32 STO+ ST X 33 FS? 48 34 X^2 35 ÷ 36 COMPLEX 37 R↓ 38 AOFF 39 END With the polynomial above: Code: 00 { 19-Byte Prgm } 01▸LBL "FX" 02 ENTER 03 ENTER 04 ENTER 05 5 06 + 07 × 08 21 09 - 10 × 11 × 12 END "FX" 4 XEQ "F'" -> 328 (-2.91e-29) 4 XEQ "F"" -> 270 (+2e-28) Cheers, Werner RE: Second derivative with complex numbers - Werner - 02-11-2021 08:53 PM I updated the routine above, because I was wrong.. Werner RE: Second derivative with complex numbers - Albert Chan - 02-12-2021 04:55 PM (02-11-2021 02:13 PM)Albert Chan Wrote:  XCas> f2(x,h) := re((f(x+h)-2*f(x)+f(x-h))/h^2) ﻿ ﻿ ﻿ ﻿ // central difference 2nd derivative ... If |h| is not too small, we avoided catastrophic cancellation, error = O(h^2) We may take advantage of it, with h = ε*√i, which make h^2 purely imaginary. With h^2 purely imaginary, we can optimized away evaluation of f(x) For h = ε*√i, all formulas below should have error = O(ε^4) f(x) ≈ re((f(x+h) + f(x-h))/2) f'(x) ≈ re((f(x+h) - f(x-h))/(2h)) f''(x) ≈ re((f(x+h) + f(x-h))/h^2) RE: Second derivative with complex numbers - peacecalc - 02-13-2021 09:00 AM hmm... a very serious contribution! [attachment=9098] RE: Second derivative with complex numbers - Werner - 02-13-2021 02:50 PM That calls for a new rewrite, simpler again. Another attempt at improvement, I take h = ε*√(2*i) = (ε,ε) - then both h and h^2 are exact (in a decimal calculator. This is the hpmuseum, after all) Code: 00 { 71-Byte Prgm } 01▸LBL 02 02 + 03 ASTO ST L 04 GTO IND ST L 05▸LBL "F'" 06 AOFF 07 GTO 00 08▸LBL "F"" 09 AON 10▸LBL 00 11 LSTO "X" 12 1ᴇ-3 13 ENTER 14 COMPLEX 15 LSTO "h" 16 XEQ 02 17 X<> "X" 18 RCL "h" 19 +/- 20 XEQ 02 21 FC? 48 22 +/- 23 RCL+ "." 24 RCL "h" 25 FC? 48 26 STO+ ST X 27 FS? 48 28 X^2 29 ÷ 30 COMPLEX 31 R↓ 32 AOFF 33 END Still, for F', I favour the f'(x) = Im(f(x+ih))/h formula, as it needs only one function evaluation. Then the routines become: Code: 00 { 77-Byte Prgm } 01▸LBL 02 02 + 03 ASTO ST L 04 GTO IND ST L 05▸LBL "F"" 06 LSTO "X" 07 1ᴇ-3 08 ENTER 09 COMPLEX 10 LSTO "h" 11 XEQ 02 12 X<> "X" 13 RCL "h" 14 +/- 15 XEQ 02 16 RCL+ "X" 17 RCL "h" 18 X^2 19 ÷ 20 COMPLEX 21 R↓ 22 RTN 23▸LBL "F'" 24 1ᴇ-99 25 COMPLEX 26 ASTO ST L 27 XEQ IND ST L 28 COMPLEX 29 1ᴇ99 30 × 31 END Cheers, Werner