# Derivative of real functions including Re and Im

When deriving functions using Re, Im or Arg (and probably some other functions as well), Mathematica gives nonsensical results, e.g.

Assuming[Element[x, Reals], D[Re[Gamma[I x]], x]]
(*
==> I Gamma[I x] PolyGamma[0, I x] Re'[Gamma[I x]]
*)

Of course the wanted result would be -Im[Gamma[I x] PolyGamma[0, I x]]. Now in this case, it is easy to resolve, but often those functions are deeply inside a complicated expression. Therefore my question is:

Is there any way to get reasonable derivatives for real functions of real parameters which involve Re, Im etc. of complex expressions depending on those variables?

-
How would you define the derivative of Im, Arg, and Re? Do these functions have a derivative? –  Searke Apr 1 '12 at 18:57
@Searke: If x is a real variable, then D[Re[f[x]],x] is Re[f'[x]], D[Im[f[x]],x] is Im[f'[x]] and D[Arg[f[x]],x] is Im[f'[x]/f[x]]. I don't think this can be handled by giving a value to Re'[f[x]], Im'[f[x]] and Arg'[f[x]] (but if I'm wrong and it can, that would allow to fix up the expression using simple replacement rules). –  celtschk Apr 1 '12 at 19:11
Leonid Shifrin´s answer below is correct. The real issue here is conceptual. If there is no definition for the derivative of of these functions, then you cannot take their derivative. –  Searke Apr 1 '12 at 19:55
@Searke: I do not want to take the derivative of Re, Im and Arg, I want to take the derivative of functions of real variables containing them. And those derivatives are defined (as long as the functions are otherwise well-behaved). –  celtschk Apr 1 '12 at 19:58
How would you handle the fact that the chain rule from calculus no longer holds? D[Im@f@x,x] = Im´[f[x]] f´[x] –  Searke Apr 1 '12 at 20:20

The functions Re and Im (just as Conjugate) don't satisfy the Cauchy-Riemann differential equations and are therefore not analytic. That means their derivative is not uniquely defined in the complex plane. That's the reason why Re' and Im' can't be simplified.

Therefore, we have to be more specific about how we want the limit to be done that corresponds to the desired derivative. The cleanest way of doing that is this simple replacement for the derivative:

Limit[(Gamma[I (x + ε)] - Gamma[I x])/ε, ε -> 0]

$i \Gamma (i x) \psi ^{(0)}(i x)$

Now I'll try to apply this to the real part instead:

Assuming[Element[x, Reals], Limit[(Re[Gamma[I (x + ε)]] - Re[Gamma[I x]])/ε, ε -> 0]]

$-\Im(\psi ^{(0)}(i x)) \Re(\Gamma (i x))-\Im(\Gamma (i x)) \Re(\psi ^{(0)}(i x))$

So we don't get any of the pesky Re' here. For those who don't like the $\LaTeX$ format, here is the real output:

-Im[PolyGamma[0, I x]] Re[Gamma[I x]] - Im[Gamma[I x]] Re[PolyGamma[0, I x]]

Edit: generalization

Motivated by the discussion of Heike's answer, I wrote a definition of a directional derivative in an arbitrary direction in the complex plane. I did this to show that such a definition can be made without in any way modifying SystemOptions (see in particular "More Information").

dirDeriv[f_, var_, angle_: 0] := Module[{g, x},
g = f /. var -> x;
Assuming[x ∈ Reals,
D[g, x] /. (h_'[p_] :> h'[p]/D[p, x]) /. (h_'[p_] :>
Limit[(h[p /. x -> x + ε Exp[I angle]] -
h[p])/ε/Exp[I angle], ε -> 0]) /.
x -> var
]
]

The first argument is the function to be differentiated, the second argument is the name of the independent variable, which will be treated as a real number. The third (optional) argument is the phase angle of the line in the complex plane along which the limit for the derivative is taken. By default this angle is zero corresponding to the real axis (the result depends on the angle only if the first argument is non-analytic).

Here is how to apply this function as a replacement for the standard derivative:

dirDeriv[Re[Gamma[I x]], x]

$-\Im(\psi ^{(0)}(i x)) \Re(\Gamma (i x))-\Im(\Gamma (i x)) \Re(\psi ^{(0)}(i x))$

A simpler example is

dirDeriv[Re[Exp[I x^2]], x]

$-2 x \sin \left(x^2\right)$

The function works by doing the standard derivative and then looking for any occurrences of Re', Im' and others (in fact, anything looking like h_'[p_]). These patterns are then replaced by first undoing the chain rule that Mathematica automatically applies (that's the division by D[p, x]), and then calculating the directional derivative as a limit analogous to the example above.

With this function one can easily explore the direction dependence of the derivative for non-analytic functions. For example,

Table[dirDeriv[Re[x], x, angle], {angle, 0, Pi, Pi/4}]

{1, 1/2 - I/2, 0, 1/2 + I/2, 1}

Edit 2: numerical functions

The above discussion concerns symbolic directional derivatives of possibly non-analytic functions. The situation is in some ways much simpler if the goal is to do numerical differentiation, i.e., calculate the derivative of a function f[x] which is numeric when the argument x is numeric.

For that case, one can go straight to the numerical approach by doing this:

Needs["NumericalCalculus`"];
ND[Re[Gamma[I x]], x, 1]

With the numerical derivative ND, you can for example make plots directly. Here is a comparison of the symbolic result with the numerical one:

Plot[{
-Im[PolyGamma[0, I x0]] Re[Gamma[I x0]] -
Im[Gamma[I x0]] Re[PolyGamma[0, I x0]],
ND[Re[Gamma[I x]], x, x0]},
{x0, 0, 5}, PlotStyle -> {Thin, Dotted}]

-
That looks like a solution! However for some reason packaging that into a function doesn't seem to work: After Dr[expr_, var_]:=Assuming[Element[var, Reals],Limit[((expr/.var->(var+eps))-expr)/eps, eps->0]] the expression Dr[Re[Gamma[I x]], x] gives 0. Any idea why? –  celtschk Apr 1 '12 at 21:19
Ok, I've found a solution: Just wrap the whole thing in With[{e = expr}, ...] and use e inside. With that it works just great! –  celtschk Apr 1 '12 at 21:26
@celtschk Hm, I tried your original function and it worked, too... maybe you had some global variable floating around. Anyway, looks like you figured out how to do the replacement. –  Jens Apr 1 '12 at 21:27
Nice solution, +1. Looks obvious in retrospect, I wish I had thought about it. I think a much simpler way (w.r.t. Cauchy - Riemann conditions) of seeing that these functions are non-analytic is to notice that they necessarily depend on both z and zbar, while CR equations are a way of saying that the function depends on z only (in addition to the existence of the derivatives). –  Leonid Shifrin Apr 1 '12 at 22:00
@LeonidShifrin Thanks for leaving something for me to do. Usually when I get to this site you and a bunch of others are already done with all the fun stuff... Whether one considers the CR equations simple is probably a matter of taste. To me, they provide the most straightforward condition. But then again, I just taught them last term, so I'm biased. –  Jens Apr 2 '12 at 3:34

Since the derivative of Re[f[x]] and Im[f[x]] with respect to a real variable x is Re[f'[x]] and Im[f'[x]], respectively, one strategy for calculating the derivatives would be exclude Re and Im when differentiating the expression and then to replace any occurrences of D[Re[a], b] and D[Im[a], b] with Re[D[a, b]] and Im[D[a, b]].

First we need to prevent Re and Im from being differentiated. We do this by adding them to the list given by SystemOptions["DifferentiationOptions" -> "ExcludedFunctions"] using SetSystemOptions:

excl = SystemOptions["DifferentiationOptions" -> "ExcludedFunctions"][[1, 2, 1, 2]];
SetSystemOptions["DifferentiationOptions" ->
"ExcludedFunctions" -> Union[Join[excl, {Re, Im}]]];

Now, when you're trying to differentiate an expression involving Re and Im Mathematica will apply the chain rule but leaves expressions of the form D[Re[f[x], x] unevaluated, e.g.

D[Sin[Re[I x^2 + 5 x]], x]

The next step is to replace the derivatives of the real or imaginary part of a function with the real or imaginary part of the derivative of the function:

deriv[exp_, x__] := D[exp, x] /. {HoldPattern[D[Re[a_], b__]] :> Re[D[a, b]],
HoldPattern[D[Im[a_], b__]] :> Im[D[a, b]]}

Example

deriv[Re[Gamma[I x]], x]

(* ==>  -Im[Gamma[I x] PolyGamma[0, I x]] *)
-
A very interesting solution. It has the advantage that it also nicely works with symbolic functions (deriv[Re[f[x]],x] gives Re[f'[x]]), however as is it doesn't work with nested Re/Im, e.g. deriv[Re[f[Re[f[x]]]],x] gives Re[D[Re[f[x]], x] f'[Re[f[x]]]]. However, replacing /. with //. seems to fix that. It has the disadvantage that you modify global state, but since I cannot imagine any case where I'd want Re' or Im', that's acceptable. Another nice feature is that it's extendable; you can also treat other (e.g. self-written) functions specially. I'll therefore move accept. –  celtschk Apr 4 '12 at 20:21
But is it really worth it having no more ability to simplify D[Re[z^2], z]? Not sure about that - but have to go now... –  Jens Apr 4 '12 at 22:27
@Heike With your modification, there is a significant change in the behavior of D[Re[Exp[I x]], x] /. x -> 1. This cannot be good. –  Jens Apr 5 '12 at 4:48
@Jens If you're worried about possible side effects you could put the SetSystemOptions statement in the definition of deriv, do the derivation, and reset the SystemOptions to their original state. –  Heike Apr 5 '12 at 9:14
@Heike That won't solve it. Try deriv[D[Re[Exp[I x z]], z] /. z -> 1, x] –  Jens Apr 5 '12 at 14:06

Ok, Jens' directional derivative functions convinced me that fixing after the fact is indeed possible; however it doesn't give nice results for non-specified symbolic functions. On the other hand, Heike's solution handles those quite nicely, but needs to modify global state, which I consider acceptable for this case (but not everyone does, as Jens' comments prove), but nevertheless solutions without modifying global state are IMHO generally preferrable. By combining elements from both approaches I've come to the following solution which seems to work quite well and also treats the cases of Abs and Arg:

deriv[expr_,var_]:=
D[expr,var]//.{Re'[e_]:>Re[D[e,var]]/D[e,var],
Im'[e_]:>Im[D[e,var]]/D[e,var],
Arg'[e_]:>Im[D[e,var]/e]/D[e,var],
Abs'[e_]:>(Re[e] Re[D[e,var]] + Im[e] Im[D[e,var]])/
(Abs[e] D[e,var])}

With this I get

deriv[Re[Gamma[I x]],x]
(*
==> -Im[Gamma[I x] PolyGamma[0, I x]]
*)
deriv[Re[f[Re[f[x]]]],x]
(*
==> Re[Re[f'[x]] f'[Re[f[x]]]]
*)
%//Simplify
(*
==> Re[f'[x]] Re[f'[Re[f[x]]]]
*)
-
Before I accept my own answer, I'll wait whether anyone finds a flaw. –  celtschk Apr 10 '12 at 11:00
OK, I think I have waited long enough :-) –  celtschk Apr 26 '12 at 16:54