Numerical differentiation for functions of complex numbers

Specific Domains Numerics
1

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.

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.

3

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.

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.

5

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.

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.

8

See Support for complex numbers · Issue #157 · JuliaDiff/ForwardDiff.jl · GitHub. 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
1 Like
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
1 Like
10

And that assumes f is analytic, correct (since it’s just the derivative along the real axis)?

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.

2 Likes
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”?

13

Yes, my brain was just being silly it seems.

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.

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.

1 Like
16

Hello, this does not seem to work for functions like sqrt. complex_diff(x -> sqrt(x), 1.0 + 0im) gives an error

StackOverflowError:

Stacktrace:
 [1] Complex at ./complex.jl:12 [inlined] (repeats 2 times)
 [2] float at ./complex.jl:1016 [inlined]
 [3] sqrt(::Complex{ForwardDiff.Dual{Nothing,Float64,1}}) at ./complex.jl:506 (repeats 79984 times)

The same error also happens if I try ForwardDiff.derivative(x -> abs(sqrt(x + 0im)), 1.0). I am curious to know how this can be done for functions involving sqrt (for a given branch). Thank you.

17

This thread is over 3 years old. Please don’t resurrect ancient threads.

A recent discussion suggested using ForwardDiff2 for differentiation of complex-valued functions.