NDSolve and differentiation of Abs

I have found several questions about the derivative of Abs and how it is not defined in the complex plane. What I have not found yet is a precise and simple workaround for when numerical calculations mix Abs and InterpolatingFunction.

Example:

r = NDSolve[{I k'[s] == 20 (s k[s] - (1 - s) (1 - k[s]) (1 + k[s])),
k[0] == 1/2}, k, {s, 0, 1}];
fun = (Abs[k[s]] /. r)[[1]];


Plotting fun'[s] is not trivial. My best solution so far:

d = 1/200;
data = ParallelTable[{ss, fun /. s -> ss}, {ss, 0, 1, d}] // N;
interpol = Interpolation[data];
Plot[interpol'[s], {s, 0, 1}, PlotRange -> All,
PlotLabel -> "Interpol"]


It works, but is quite ugly though, evaluating one InterpolatingFunction in order to create another. Also, the smaller the value of d, the more precise and slower the execution...

Would anyone point me to a more elegant and efficient way to do the same task?

Final note: The cake goes for MarcoB! J.M.'s method was more elegant (maybe even faster), but Marcos' can work up to any order with simple modifications. Thank you both for your excellent contributions.

You could try using numerical differentiation from the Numerical Calculus package:

Clear[r, fun]
Needs["NumericalCalculus"]

r = NDSolve[
{I k'[s] == 20 (s k[s] - (1 - s) (1 - k[s]) (1 + k[s])), k[0] == 1/2},
k, {s, 0, 1}
];
fun[s_] := Abs[(k /. r[[1]])[s]];

Plot[
{10 fun[t], ND[fun[s], s, t, Terms -> 20]},
{t, 0, 1},
PlotLegends -> {"scaled function", "numerical derivative"}
]


Note that I scaled up the values of the function 10x for convenience, so it would be visible on the same range as the derivative. The choice of the number of terms to include in the ND expression was somewhat arbitrary; the default of 7 terms, however, was rather too low for this function and led to artifacts.

• Excellent!!! Didnt know about the Terms option and was just about to ask. – ivbc Jun 16 '16 at 20:54
• @ivbc Glad it helped! – MarcoB Jun 16 '16 at 21:13
• I just noticed that we defined fun in somewhat different ways. Is that just style or something deeper? – ivbc Jun 17 '16 at 13:16
• @ivbc It was mostly a remnant of other trials I had made before when I was toying with the Terms options. I had thought that perhaps explicitly exposing Abs to ND might help if ND carries out any symbolic pre-evaluation (like NIntegrate would). It turns out that it doesn't matter, but I forgot to revert the changes when I posted the answer. I did so now. Defining fun as a SetDelayed (:=) expression protects me in case I found a need to change the name of the independent variable down the road. – MarcoB Jun 17 '16 at 16:27

Alternatively, you could do a little complex number algebra for the purpose:

kf = NDSolveValue[{I k'[s] == 20 (s k[s] - (1 - s) (1 - k[s]) (1 + k[s])),
k[0] == 1/2}, k, {s, 0, 1}];

Plot[{10 Abs[kf[t]], Re[kf[t] Conjugate[kf'[t]]]/Abs[kf[t]]}, {t, 0, 1},
PlotLegends -> {"scaled function", "numerical derivative"}]


• Quite elegant trick for skipping the derivative of the Abs!!! Why did you use NDSolveValue instead of NDSolve? – ivbc Jun 17 '16 at 18:24
• @ivbc, yes, I could have used NDSolve[], but I didn't wan't to have to bother with an additional use of First[] and /.. – J. M. will be back soon Jun 17 '16 at 22:54

I have found yet another way to solve this. Based on J.M.'s answer and another post on defining the derivative of Abs for Real valued functions (will look it up later to add here).

The following code works:

Derivative[1][Abs][x_] := Re[  Conjugate[x] D[x, s]  ]/Abs[x]/D[x, s];
Plot[D[Abs[k[s]], s] /. k -> kf  /. s -> ss, {ss, 0, 1}]


It is worth noticing that now we can easily plot derivatives of compound functions of Abs.

Plot[D[1/(1 + Abs[k[s]]), s] /. k -> kf  /. s -> ss, {ss, 0, 1}]


We can go a bit further and add the second derivative as well!

Derivative[2][Abs][
x_] := (Re[Conjugate[x] D[x, {s, 2}]]/Abs[x] +
Im[Conjugate[x] D[x, s]]^2/(Abs[x]^3))/D[x, {s, 2}];
Derivative[1][Abs][x_] := Re[Conjugate[x] D[x, s]]/Abs[x]/D[x, s];
Derivative[1][Re][x_] := Re[D[x, s]]/D[x, s];
Derivative[1][Conjugate][x_] := Conjugate[D[x, s]]/D[x, s];
Plot[D[Abs[k[s]],{s,2}] /. k -> kf  /. s -> ss, {ss, 0, 1}]


There is one problem with this solution though: D[x, s] means that we can only derivate Abs in respect to s. I dont know how to generalize this and any help is welcome.

Also, in my formula for the first derivative of Abs, there is one extra derivative dividing the whole expression. It is different from the expression proposed by J.M. only because of the chain rule, because Mathematica needs the definition of dAbs[k[s]]/dk, and not dAbs[k[s]]/ds.

Note:All mathematical formulas presented here can be verified by using the complex numbers notation k=a*E^(Ib).

• To assist you with the derivation: ponder on the result of D[Sqrt[r[t]^2 + i[t]^2], t]`. – J. M. will be back soon Jun 18 '16 at 14:07