(Insanely-) Higher Order Derivatives

Specific Domains Numerics
21

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.

2 Likes
22

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.

23

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)

24

Could you clarify this formula? What is w_i? What is being plugged into the trapezoidal rule?

25

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.

1 Like
26

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.

1 Like
27

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).

28

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?
29

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 on the trapezoidal rule for contour integrals and similar cases, for example.)

1 Like
30

I have the same impression.

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

2 Likes
31

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.

6 Likes
32

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

The quick-and-dirty script I used was

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.

5 Likes
33

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

5 Likes
34

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

35

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 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…)

2 Likes
36

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).

1 Like
37

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.

1 Like