# (Insanely-) Higher Order Derivatives

**URL:** <https://discourse.julialang.org/t/insanely-higher-order-derivatives/104455>\
**Category:** Numerics\
**Tags:** question\
**Created:** [October 1, 2023, 8:12am UTC](https://discourse.julialang.org/t/insanely-higher-order-derivatives/104455 "2023-10-01T08:12:55Z")\
**Posts on this page:** 17\
**Page:** 2

<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:** [October 5, 2023, 9:23pm UTC](https://discourse.julialang.org/t/insanely-higher-order-derivatives/104455/21 "2023-10-05T21:23:40Z")

</div>

> [@dlfivefifty](#):
>
> The issue though is that I believe Bornemann chooses the contour by optimising a condition number coming from the trapezium rule. I don’t know if you can do this with QuadGK.

From the paper, the condition number is given in terms of the mean of the absolute value around the contour, which you can just as easily compute with QuadGK or any integration method (especially since you probably only need the condition number to low accuracy).

On the other hand, if you adaptively choose the contour, you probably end up with a contour that is not close to a pole(?), in which case the trapezoidal rule should be fine.

---

<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:** [October 5, 2023, 9:29pm UTC](https://discourse.julialang.org/t/insanely-higher-order-derivatives/104455/22 "2023-10-05T21:29:59Z")

</div>

I forget the details but if I recall correctly most functions are “exponential like”: actually only a small part of the contour contributes to the final result. So its possible you are right that QuadGK still uses less function evaluations.

---

<div class="post-metadata">

**Author:** ![Veenty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/veenty/32/50940_2.png) [@Veenty](https://discourse.julialang.org/u/Veenty)\
**Post date:** [October 5, 2023, 9:53pm UTC](https://discourse.julialang.org/t/insanely-higher-order-derivatives/104455/23 "2023-10-05T21:53:36Z")

</div>

To be fair, there are a couple of nice tricks that can be done to compute that integral still using trapezoidal rule.

A good one is noting that if f is analytic, then

f(z) ~ \frac{\sum\_{i} \frac{f(w\_i)}{w\_i - z}}{\sum\_{i} \frac{1}{w\_i - z} }

And the convergence of this is magically better for close evaluations (d(z, \Gamma) \<\< 1)

---

<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:** [October 6, 2023, 12:23am UTC](https://discourse.julialang.org/t/insanely-higher-order-derivatives/104455/24 "2023-10-06T00:23:29Z")

</div>

> [@Veenty](#):
>
> A good one is noting that …

Could you clarify this formula? What is w\_i? What is being plugged into the trapezoidal rule?

---

<div class="post-metadata">

**Author:** ![Veenty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/veenty/32/50940_2.png) [@Veenty](https://discourse.julialang.org/u/Veenty)\
**Post date:** [October 6, 2023, 12:44am UTC](https://discourse.julialang.org/t/insanely-higher-order-derivatives/104455/25 "2023-10-06T00:44:09Z")

</div>

Let \gamma:[0,1] \to \mathbb{C} be a smooth parametrization of the curve \Gamma. Let w\_i = \gamma(\frac{i}{N}) with i = 0... N-1. From the trapezoidal approximation and cauchy formula it follows that

f(z) = \frac{1}{2 \pi i }\int\_{\Gamma} \frac{f(w)}{w - z} dw = \frac{\int\_{\Gamma} \frac{f(w)}{w - z} dw}{\int\_{\Gamma} \frac{1}{w - z} dw} \approx \frac{ \frac{1}{N}\sum\_{i=0}^{N-1} \frac{f(w\_i )}{w\_i - z} \gamma^\prime(\frac{i}{N})}{\frac{1}{N}\sum\_{i=0}^{N-1} \frac{1}{w\_i - z} \gamma^\prime(\frac{i}{N}) }

What its almost magic about this representation is that the errors for z close the the boundary magically cancel out, and this expression ends up being an extremely good approximation.

---

<div class="post-metadata">

**Author:** ![Veenty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/veenty/32/50940_2.png) [@Veenty](https://discourse.julialang.org/u/Veenty)\
**Post date:** [October 6, 2023, 1:12am UTC](https://discourse.julialang.org/t/insanely-higher-order-derivatives/104455/26 "2023-10-06T01:12:35Z")

</div>

But the question still remains, what are you after that requieres so many derivatives? Is it a complicated function that the taylor expansion is easier to evaluate? Has the derivative itself has any meaning? Either way, taking 4 derivatives is already super ill-condition problem numerically, so even if your function is analytic I would be shocked if you are able to get 1 decimal correctly if you are taking 1000 derivatives.

---

<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:** [October 6, 2023, 1:24am UTC](https://discourse.julialang.org/t/insanely-higher-order-derivatives/104455/27 "2023-10-06T01:24:49Z")

</div>

> [@Veenty](#):
>
> What its almost magic about this representation is that the errors for z close the the boundary magically cancel out, and this expression ends up being an extremely good approximation.

I don’t see how it helps with the example I gave, assuming you _don’t_ know the location of the pole. If you know the location z of the pole(s), I agree that there are lots of useful tricks (both for the trapezoidal rule and for other numerical integration schemes).

---

<div class="post-metadata">

**Author:** ![Syx\_Pek](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/syx_pek/32/6364_2.png) [@Syx\_Pek](https://discourse.julialang.org/u/Syx_Pek)\
**Post date:** [October 6, 2023, 6:13am UTC](https://discourse.julialang.org/t/insanely-higher-order-derivatives/104455/28 "2023-10-06T06:13:47Z")

</div>

Perhaps I’m being dense here, but

1. Why is taking this summation any better than dividing than 2\pi i? The denominator here is effectively 2\pi i. Is it just that \frac 1N cancels?
2. Perhaps this is a restatement of Steven’s concern, the choice of \Gamma needs to ensure that f is holomorphic (equivalently, analytic) in the interior of \Gamma. Do you have a general scheme to choose \Gamma?

---

<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:** [October 6, 2023, 10:02am UTC](https://discourse.julialang.org/t/insanely-higher-order-derivatives/104455/29 "2023-10-06T10:02:43Z")

</div>

> [@Syx\_Pek](#):
>
> Why is taking this summation any better than dividing than 2\pi i? The denominator here is effectively 2\pi i.

It is only 2\pi i in the limit N \to \infty. For finite N there are errors, and these errors can grow large when the (known) pole z is close to the contour \Gamma; from what @Veenty said, these finite-N errors somehow mostly cancel between the numerator and denominator.

(It would be nice to have a reference; this technique isn’t mentioned in [Trefethen’s paper](https://people.maths.ox.ac.uk/trefethen/sirev56-3_385.pdf) on the trapezoidal rule for contour integrals and similar cases, for example.)

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [October 6, 2023, 11:32am UTC](https://discourse.julialang.org/t/insanely-higher-order-derivatives/104455/30 "2023-10-06T11:32:28Z")

</div>

> [@Veenty](#):
>
> But the question still remains, what are you after that requieres so many derivatives?

I have the same impression.

> [@NitinGupta](#):
>
> two words : “Theoretical Physics” 😂

is a pithy answer that does not tell us anything. Knowing the application would make it easier to help.

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [October 6, 2023, 12:23pm UTC](https://discourse.julialang.org/t/insanely-higher-order-derivatives/104455/31 "2023-10-06T12:23:46Z")

</div>

I could imagine that such a pithy reply might be okay-ish if people in the field knew what he was talking about, but I strongly suspect that if you did a survey of theoretical physicists asking them if they knew of a single use for the 1000th order derivative of a function, the number who answered yes would be around the same number as those who believed the world was flat.

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [October 6, 2023, 1:26pm UTC](https://discourse.julialang.org/t/insanely-higher-order-derivatives/104455/32 "2023-10-06T13:26:56Z")

</div>

I think there are still a couple of errors in this. You wrote dz instead of dw and the dw isn’t included in the sums. But I tried what I think you intended to evaluate \cos(z) at the point where @stevengj placed the pole in his example using the unit circle as the contour. And I have to agree with you that this does seem quite magical. The plot I got was  
 ![trapezoidal](https://global.discourse-cdn.com/julialang/original/3X/0/1/01e0cb1e79b2a7a9a210b2ba83917ebfc010f95b.png)

The quick-and-dirty script I used was

```julia
using Plots

function trap(f, z, n)
  sum = 0.0im
  for j in 0:(n - 1)
    w = exp(im * 2 * pi * j / n)
    sum += f(w) * w / (w - z)
  end
  return sum / n
end

function trap_ratio(f, z, n)
  sum_num = 0.0im
  sum_den = 0.0im
  for j in 0:(n - 1)
    w = exp(im * 2 * pi * j / n)
    sum_num += f(w) * w / (w - z)
    sum_den += w / (w - z)
  end
  return sum_num / sum_den
end

z = 0.9790 + 0.198im
f(z) = cos(z)

ns = [2:2:50; (k -> round(10^(k/4))).(8:20)]
results = (n -> abs(trap(f,z,n)-f(z))).(ns)
results_ratio = (n -> abs(trap_ratio(f,z,n)-f(z))).(ns)

plot(
  ns,
  results,
  yaxis = :log10,
  xaxis = :log10,
  yticks = (k -> 10.0^(-2 * k)).(-1:16),
  marker = :circle,
  markersize = 2,
  label = "Trapezoidal",
  xlabel = "points",
  ylabel = "error",
)
plot!(ns, results_ratio, marker=:diamond, markersize=2, label="Trapezoidal Ratio")

```

I’d like to second the request for a reference if you have it. I’d really like to understand how this works.

---

<div class="post-metadata">

**Author:** ![Veenty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/veenty/32/50940_2.png) [@Veenty](https://discourse.julialang.org/u/Veenty)\
**Post date:** [October 6, 2023, 3:47pm UTC](https://discourse.julialang.org/t/insanely-higher-order-derivatives/104455/33 "2023-10-06T15:47:25Z")

</div>

Its on page 39 of the PDF of that paper! Trefethen has an explanation in the circle using rational functions.  
The original paper is this [one](https://link.springer.com/article/10.1007/BF01931287)

---

<div class="post-metadata">

**Author:** ![Veenty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/veenty/32/50940_2.png) [@Veenty](https://discourse.julialang.org/u/Veenty)\
**Post date:** [October 6, 2023, 4:00pm UTC](https://discourse.julialang.org/t/insanely-higher-order-derivatives/104455/34 "2023-10-06T16:00:10Z")

</div>

My mistake, I corrected it to include the Jacobian now!

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [October 6, 2023, 5:53pm UTC](https://discourse.julialang.org/t/insanely-higher-order-derivatives/104455/35 "2023-10-06T17:53:04Z")

</div>

Interesting. So it has a straight-forward (if perhaps obscure) explanation based on the fact that it’s really barycentric Lagrange interpolation. I skimmed over that page before responding previously, but didn’t notice what the Matlab code was actually doing. (Berrut and Trefethen’s paper on barycentric interpolation [https://people.maths.ox.ac.uk/trefethen/barycentric.pdf](https://people.maths.ox.ac.uk/trefethen/barycentric.pdf) is also a good read.)

Unfortunately it looks like if you want f^{(k)}(z) for k\geq 1 and have the factor 1/(w-z)^{k+1} in the integral, you will still have the delay in exponential convergence when z is close to the contour. (Unless there is some further trick…)

---

<div class="post-metadata">

**Author:** ![Veenty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/veenty/32/50940_2.png) [@Veenty](https://discourse.julialang.org/u/Veenty)\
**Post date:** [October 6, 2023, 9:56pm UTC](https://discourse.julialang.org/t/insanely-higher-order-derivatives/104455/36 "2023-10-06T21:56:56Z")

</div>

In the original paper (that I also linked) have general formulas for higher order derivatives. However, it seems that you loose one digit of precision for very derivative that you take (which is quite standard with derivatives).

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [October 6, 2023, 10:29pm UTC](https://discourse.julialang.org/t/insanely-higher-order-derivatives/104455/37 "2023-10-06T22:29:58Z")

</div>

I was hoping for the further trick. I guess I should look at that paper. Unfortunately I don’t have institutional online access. I’ll have to put a little more work into getting a copy, but it sounds worthwhile.

[Previous page](https://discourse.julialang.org/t/insanely-higher-order-derivatives/104455.md?page=1)