Does anyone know of an available package which supports differentiating complex numbers? This has come up when trying to make stiff solvers for ODEs compatible with complex numbers since numerical Jacobians usually need to be calculated (unless the user provides a function for it), but when checking the standard solutions (Calculus.jl, ForwardDiff.jl) I came up empty handed.

# Numerical differentiation for functions of complex numbers

**tshort**#2

You could wrap zvode. Itâs one Fortran subroutine, so it should be easy to wrap.

As a long shot, you might also look into the ApproxFun package. In the Matlab world, the analogous package is chebfun, and I think itâs been used for ODEs with complex numbers. See the example here.

Thank you for the response, I donât think either of these provide what I am looking for.

Iâd rather try to extend the native Julia routines to complex numbers. Weâre already âfaster than Fortranâ in many cases and just need to take complex derivatives in order to get into that domain, so I donât think thereâs a reason to give up so soon.

Does ApproxFun handle arbitrary functions? Every example I see is suspiciously simple and on mathematical functions. Will it differentiate things with for loops and conditionals in it? I would like to start using it more, but alas without docs and without prior knowledge of chebfun I am scared. But if someone would guide me into that strange new world I would definitely like to start using it in the DiffEq universe.

**lbenet**#4

TaylorSeries.jl does this to some extent:

```
julia> using TaylorSeries
julia> x,y = set_variables(Complex128, "x y", order=5)
2-element Array{TaylorSeries.TaylorN{Complex{Float64}},1}:
( 1.0 ) x + đȘ(âxââ¶)
( 1.0 ) y + đȘ(âxââ¶)
julia> f(x,y) = [x^3-2*im*y^2+x-y, x-x*y-y+1im]
f (generic function with 1 method)
julia> jacobian(f(x,y))
2Ă2 Array{Complex{Float64},2}:
1.0+0.0im -1.0+0.0im
1.0+0.0im -1.0+0.0im
julia> jacobian(f(x,y), [complex(2,1), complex(-1,-3)])
2Ă2 Array{Complex{Float64},2}:
10.0+12.0im -13.0+4.0im
2.0+3.0im -3.0-1.0im
```

I hope this helps.

But this too has a strong restriction on the types of functions it can differentiate, right? These tools are nice and I will like to incorporate them more, but I am looking for something where I can give it a very general class of Julia-defined functions and get a derivative out for complex numbers, like what ForwardDiff.jl and Calculus.jl offer for real-valued functions.

**lbenet**#6

You are right, TaylorSeries.jl doesnât support all functions implemented in Julia, yet. Simply put, we had no motivation to implement them. You can certainly implement what ever you need and TaylorSeries lacks, and send a PR.

**jrevels**#8

See https://github.com/JuliaDiff/ForwardDiff.jl/issues/157. ForwardDiffâs `Dual`

type is capable of complex differentiation, but I havenât thought too deeply about the API for doing so.

Hereâs an example if you want to play around with it:

```
julia> using ForwardDiff: Dual, partials
# a and b are my real and imag components, respectively
julia> a, b = rand(2)
2-element Array{Float64,1}:
0.125707
0.388673
# let's get the derivative of sin w.r.t. a
# (and, by way of Cauchy-Riemann, w.r.t. b - thanks @dpsanders)
julia> d = sin(Dual(a, one(a)) + Dual(b, zero(b))*im)
Dual(0.13496604124570097,1.0679948433847894) + Dual(0.39538861882848864,-0.0499665676921812)*im
julia> cos(a + b*im) == partials(real(d), 1) + partials(imag(d), 1) * im
true
```

**dpsanders**#9

You can wrap this up as

```
function complex_diff(f, z)
a, b = reim(z)
z_dual = Dual(a, one(a)) + im * Dual(b, zero(b))
deriv = f(z_dual)
return partials(real(deriv), 1) + im * partials(imag(deriv), 1)
end
```

**ChrisRackauckas**#10

And that assumes `f`

is analytic, correct (since itâs just the derivative along the real axis)?

**stevengj**#11

ApproxFun works on arbitrary functions `f(x)`

, and doesnât care how you compute `f`

: it just uses the function as a âblack boxâ that produces a value given `x`

. It works best with smooth functions, but if you want to compute derivatives then I assume that your function is at least differentiable.

**dpsanders**#12

Doesnât the definition of a complex function being differentiable at a point include implicitly that the derivative is the same in each âdirectionâ?

**ChrisRackauckas**#14

Thanks for the clarification. Iâll have to take another look at the ApproxFun then. I have a lot of questions about it but will put it in a different thread.

**dlfivefifty**#15

This bit of code also calculates the Jacobian as in @lbenet example:

```
using ApproxFun
d=Interval(0,2+im)*Interval(-1-3im,0)
x,y=Fun(d)
f=[x^3-2im*y^2+x-y;x-x*y-y+1im]
Dx=Derivative(d,[1,0])+0im # + 0im is workaround for a bug
Dy=Derivative(d,[0,1])+0im
J=[Dx*f[1] Dy*f[1];
Dx*f[2] Dy*f[2]]
map(g->g(2+im,-1-3im),J)
```

Itâs not quite as elegant as it should be yet, though.