# Numerical derivative of numerical function

I am trying to calculate the derivative of a numerically defined Hamiltonian. This is a minimal working example (on a simple example):

H[kx_?NumericQ,
ky_?NumericQ] := {{0, 0 (I kx + ky) , 0, 0}, {(-I kx + ky) , 0, 0,
0}, {0, 0, 0, I ((kx - I ky) )}, {0, 0, I ((kx + I ky) ), 0}}

v1[kx_?NumericQ, ky_?NumericQ] :=
v1[kx, ky] = Eigensystem[H[kx, ky]][[2, 1]];
v1c[kx_?NumericQ, ky_?NumericQ] :=
v1c[kx, ky] = Refine[Conjugate[Eigensystem[H[kx, ky]]][[2, 1]]];

v1n[kx_?NumericQ, ky_?NumericQ] :=
N[v1[kx, ky]/Sqrt[v1c[kx, ky] . v1[kx, ky]]];
v1cn[kx_?NumericQ, ky_?NumericQ] :=
v1c[kx, ky]/Sqrt[v1c[kx, ky] . v1[kx, ky]];

dumv1n[u_?NumericQ, w_?NumericQ] :=
N[v1n[kx, ky] /. {kx -> u, ky -> w}];
Dumdv1nN[p_?NumericQ, q_?NumericQ] :=
N[D[dumv1n[u, w], u] /. {u -> p, w -> q}];
dv1n[p_?NumericQ, q_?NumericQ] := N[Dumdv1nN[p, q]];


But it doesn't calculate the derivative when I run dv1n[10,10].

• Suppose $H(t)$ happens to be a square matrix depending on $t$ with distinct eigenvalues and a unique eigenvalue of largest absolute value. Let $v(t)$ be a normalized eigenvector corresponding to the eigenvalue. Then $\alpha(t)\,v(t)$ is also a normalized eigenvector, where $\alpha$ is an arbitrary function with $|\alpha(t)|=1$. Now $\alpha(t)\,v(t)$ is not guaranteed to be differentiable or even continuous. Even if, say, both $\alpha(t)$ and $v(t)$ are differentiable, there are many different derivatives. Which of these $\alpha(t)\,v(t)$, if any, does Eigensystem compute? Commented Jun 13 at 17:52

Maybe you can use the function ND in the package "NumericalCalculus". Just modify your code of function Dumdv1nN:

Needs["NumericalCalculus"]

H[kx_?NumericQ,
ky_?NumericQ] := {{0, 0  (I  kx + ky), 0, 0}, {(-I  kx + ky), 0, 0,
0}, {0, 0, 0, I  ((kx - I  ky))}, {0, 0, I  ((kx + I  ky)), 0}}

v1[kx_?NumericQ, ky_?NumericQ] :=
v1[kx, ky] = Eigensystem[H[kx, ky]][[2, 1]];
v1c[kx_?NumericQ, ky_?NumericQ] :=
v1c[kx, ky] = Refine[Conjugate[Eigensystem[H[kx, ky]]][[2, 1]]];

v1n[kx_?NumericQ, ky_?NumericQ] :=
N[v1[kx, ky]/Sqrt[v1c[kx, ky] . v1[kx, ky]]];
v1cn[kx_?NumericQ, ky_?NumericQ] :=
v1c[kx, ky]/Sqrt[v1c[kx, ky] . v1[kx, ky]];

dumv1n[u_?NumericQ, w_?NumericQ] :=
N[v1n[kx, ky] /. {kx -> u, ky -> w}];
Dumdv1nN[p_?NumericQ, q_?NumericQ] := ND[dumv1n[u, q], u, p];
dv1n[p_?NumericQ, q_?NumericQ] := N[Dumdv1nN[p, q]];



Now it works:

dv1n[10, 10]


{0. + 0. I, 0. + 0. I, -9.99684*10^-15 + 0. I, 218.316 + 218.026 I}

• Don't be bullied by that user about the log Pearson Type III distribution. Many of us have been on the receiving end of such comments (sometimes right and sometimes - as in this case - wrong).
– JimB
Commented Jun 14 at 17:34
• But I think you will need to modify one thing. It is the log of the Log Pearson Type III has a shifted gamma distribution (analogous to the log of a LogNormal is a normal distribution.
– JimB
Commented Jun 14 at 17:42

dumv1n` accepts only numerical arguments, but you try to feed it symbols (when you attempt symbolic differentiation).