# Numerical differentiation for functions of complex numbers

**URL:** <https://discourse.julialang.org/t/numerical-differentiation-for-functions-of-complex-numbers/341>\
**Category:** Numerics\
**Tags:** complex-numbers\
**Created:** [November 16, 2016, 2:56am UTC](https://discourse.julialang.org/t/numerical-differentiation-for-functions-of-complex-numbers/341 "2016-11-16T02:56:52Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 16, 2016, 2:56am UTC](https://discourse.julialang.org/t/numerical-differentiation-for-functions-of-complex-numbers/341/1 "2016-11-16T02:56:52Z")

</div>

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](https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/110) 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.

---

<div class="post-metadata">

**Author:** ![tshort](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tshort/32/43_2.png) [@tshort](https://discourse.julialang.org/u/tshort)\
**Post date:** [November 17, 2016, 12:15pm UTC](https://discourse.julialang.org/t/numerical-differentiation-for-functions-of-complex-numbers/341/2 "2016-11-17T12:15:14Z")

</div>

You could wrap [zvode](http://netlib.sandia.gov/ode/zvode.f). It’s one Fortran subroutine, so it should be easy to wrap.

As a long shot, you might also look into the [ApproxFun](https://github.com/ApproxFun/ApproxFun.jl) 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](http://www.chebfun.org/examples/ode-linear/DynamicalSystems.html).

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 17, 2016, 12:30pm UTC](https://discourse.julialang.org/t/numerical-differentiation-for-functions-of-complex-numbers/341/3 "2016-11-17T12:30:33Z")

</div>

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

> [@tshort.rlists](#):
>
> You could wrap zvode. It’s one Fortran subroutine, so it should be easy to wrap.

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.

> [@tshort.rlists](#):
>
> 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.

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.

---

<div class="post-metadata">

**Author:** ![lbenet](https://avatars.discourse-cdn.com/v4/letter/l/35a633/32.png) [@lbenet](https://discourse.julialang.org/u/lbenet)\
**Post date:** [November 17, 2016, 1:49pm UTC](https://discourse.julialang.org/t/numerical-differentiation-for-functions-of-complex-numbers/341/4 "2016-11-17T13:49:53Z")

</div>

[TaylorSeries.jl](https://github.com/JuliaDiff/TaylorSeries.jl) does this to some extent:

```julia
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 17, 2016, 1:55pm UTC](https://discourse.julialang.org/t/numerical-differentiation-for-functions-of-complex-numbers/341/5 "2016-11-17T13:55:01Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![lbenet](https://avatars.discourse-cdn.com/v4/letter/l/35a633/32.png) [@lbenet](https://discourse.julialang.org/u/lbenet)\
**Post date:** [November 17, 2016, 2:43pm UTC](https://discourse.julialang.org/t/numerical-differentiation-for-functions-of-complex-numbers/341/6 "2016-11-17T14:43:42Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![jrevels](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jrevels/32/10393_2.png) [@jrevels](https://discourse.julialang.org/u/jrevels)\
**Post date:** [November 17, 2016, 3:43pm UTC](https://discourse.julialang.org/t/numerical-differentiation-for-functions-of-complex-numbers/341/8 "2016-11-17T15:43:58Z")

</div>

See [Support for complex numbers · Issue #157 · JuliaDiff/ForwardDiff.jl · GitHub](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
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

```

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [November 17, 2016, 6:37pm UTC](https://discourse.julialang.org/t/numerical-differentiation-for-functions-of-complex-numbers/341/9 "2016-11-17T18:37:02Z")

</div>

You can wrap this up as

```julia
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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 17, 2016, 7:55pm UTC](https://discourse.julialang.org/t/numerical-differentiation-for-functions-of-complex-numbers/341/10 "2016-11-17T19:55:21Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 17, 2016, 7:56pm UTC](https://discourse.julialang.org/t/numerical-differentiation-for-functions-of-complex-numbers/341/11 "2016-11-17T19:56:44Z")

</div>

> [@ChrisRackauckas](#):
>
> 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?

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.

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [November 17, 2016, 8:12pm UTC](https://discourse.julialang.org/t/numerical-differentiation-for-functions-of-complex-numbers/341/12 "2016-11-17T20:12:56Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 17, 2016, 8:24pm UTC](https://discourse.julialang.org/t/numerical-differentiation-for-functions-of-complex-numbers/341/13 "2016-11-17T20:24:57Z")

</div>

Yes, my brain was just being silly it seems.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 17, 2016, 8:28pm UTC](https://discourse.julialang.org/t/numerical-differentiation-for-functions-of-complex-numbers/341/14 "2016-11-17T20:28:53Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [November 19, 2016, 7:50am UTC](https://discourse.julialang.org/t/numerical-differentiation-for-functions-of-complex-numbers/341/15 "2016-11-19T07:50:03Z")

</div>

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![zwwi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zwwi/32/11932_2.png) [@zwwi](https://discourse.julialang.org/u/zwwi)\
**Post date:** [April 2, 2020, 3:46am UTC](https://discourse.julialang.org/t/numerical-differentiation-for-functions-of-complex-numbers/341/16 "2020-04-02T03:46:50Z")

</div>

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [April 2, 2020, 4:09am UTC](https://discourse.julialang.org/t/numerical-differentiation-for-functions-of-complex-numbers/341/17 "2020-04-02T04:09:55Z")

</div>

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.