# Best method for calculating the second derivative of 1D array

**URL:** <https://discourse.julialang.org/t/best-method-for-calculating-the-second-derivative-of-1d-array/69353>\
**Category:** Performance\
**Tags:** question, package, performance, differentiation, splines\
**Created:** [October 7, 2021, 8:37am UTC](https://discourse.julialang.org/t/best-method-for-calculating-the-second-derivative-of-1d-array/69353 "2021-10-07T08:37:26Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![AngusRussell](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/angusrussell/32/22226_2.png) [@AngusRussell](https://discourse.julialang.org/u/AngusRussell)\
**Post date:** [October 7, 2021, 8:37am UTC](https://discourse.julialang.org/t/best-method-for-calculating-the-second-derivative-of-1d-array/69353/1 "2021-10-07T08:37:26Z")

</div>

![Eigen_energies](https://global.discourse-cdn.com/julialang/original/3X/0/9/09abfd3edaaa2b4ce95040eaa70e8b7cffb0f88e.jpeg)

Essentially I have 1D arrays of eigenvalues as shown in the above plot and I want to find the 2nd derivative of each energy band with respect to ϵ.

My current method uses [Dierckx.jl](https://github.com/kbarbary/Dierckx.jl).

In my code, the eigenvalues are stored in a 3D array with the structure  
[_Energy band_, _value of ϵ_, _value of B_], at the moment that’s a [5, 520, 520]. Hence the loops in the code below:

```julia
function QuantumCapacitance(E, ϵ)
    Q_C = Array{Float64, 3}(undef, 5, size(E,2), size(E,3)) # Define array
    for i = 1:size(E,2) # Magetic field loop
        for j = 1:5 # Eigenvalue Loop
            itp = Spline1D(ϵ, E[j,:,i], k=2) # Value for 2nd derivative 
            for k = 1:size(E,2) # detuning, ϵ, loop 
            Q_C[j, k, i] = derivative(itp,ϵ[k],nu=2) 
            end
        end
    end
    return Q_C

```

This method however takes up a lot of time in my code. Is there a more efficient way to do this, are there any finite difference method packages that would work better here, would [Interpolations.jl](https://juliamath.github.io/Interpolations.jl/latest/) be quicker.

My reason for not using that package was earlier in my code I had to perform a cubic spline and Dierckx was much more accurate.

Thank you 🙂

---

<div class="post-metadata">

**Author:** ![cmarcotte](https://avatars.discourse-cdn.com/v4/letter/c/a3d4f5/32.png) [@cmarcotte](https://discourse.julialang.org/u/cmarcotte)\
**Post date:** [October 7, 2021, 9:35am UTC](https://discourse.julialang.org/t/best-method-for-calculating-the-second-derivative-of-1d-array/69353/2 "2021-10-07T09:35:05Z")

</div>

Interpolations.jl includes the `gradient` function, which may work (not all interpolants supplied by that package support gradient computation) and may be faster.  
However, I notice your `k` loop over the values of `ϵ`; Dierckx.jl will take a vector of values to compute the derivative at, so if you do `derivative(itp, ϵ, nu=2)` it may be faster and then you can assign with `Q_c[j,:,i] .=` to avoid allocations (I think? since the derivative will return an array).

---

<div class="post-metadata">

**Author:** ![Bardo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bardo/32/21601_2.png) [@Bardo](https://discourse.julialang.org/u/Bardo)\
**Post date:** [October 7, 2021, 8:41pm UTC](https://discourse.julialang.org/t/best-method-for-calculating-the-second-derivative-of-1d-array/69353/3 "2021-10-07T20:41:18Z")

</div>

Cubic spline gives second order accuracy, after 2nd derivation, not much is left…  
Your functions are smooth, so you can use higher order interpolants.  
But raising the order usually increases the effort and often runs into numerical issues.

One method without such issues is _barycentric interpolation_.  
There are a couple of Julia packages available, not all with working derivatives though.

---

<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 7, 2021, 9:23pm UTC](https://discourse.julialang.org/t/best-method-for-calculating-the-second-derivative-of-1d-array/69353/4 "2021-10-07T21:23:45Z")

</div>

> [@Bardo](#):
>
> One method without such issues is _barycentric interpolation_ .

You shouldn’t use high-order polynomial interpolation (barycentric or not) from equally spaced points, because you can run into [Runge phenomena](https://en.wikipedia.org/wiki/Runge%27s_phenomenon).

(Equally spaced points are fine if you do a least-square fit to a much lower-degree polynomial.)

---

<div class="post-metadata">

**Author:** ![Bardo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bardo/32/21601_2.png) [@Bardo](https://discourse.julialang.org/u/Bardo)\
**Post date:** [October 7, 2021, 10:06pm UTC](https://discourse.julialang.org/t/best-method-for-calculating-the-second-derivative-of-1d-array/69353/5 "2021-10-07T22:06:13Z")

</div>

True for the basic barycentric method [Berrut](https://people.maths.ox.ac.uk/trefethen/barycentric.pdf) et al.  
Two variants are addressing these issues:

- Barycentric interpolation via AAA
- Rational barycentric interpolation,

both in [GitHub - macd/BaryRational.jl: Barycentric rational approximation and interpolation in one dimension.](https://github.com/macd/BaryRational.jl)

> NOTE: Currently, the derivatives are not working.

And right, a [local fit](https://en.wikipedia.org/wiki/Local_regression) should be tried too, YMMV.

---

<div class="post-metadata">

**Author:** ![fedoroff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fedoroff/32/53209_2.png) [@fedoroff](https://discourse.julialang.org/u/fedoroff)\
**Post date:** [October 8, 2021, 6:37am UTC](https://discourse.julialang.org/t/best-method-for-calculating-the-second-derivative-of-1d-array/69353/6 "2021-10-08T06:37:06Z")

</div>

What about simple [finite differences](https://en.wikipedia.org/wiki/Finite_difference)?

```julia
function derivative2(x, y, i)
    Nx = length(x)
    @assert Nx >= 5
    dx = x[2] - x[1]
    if i == 1
        dy = 35*y[i] - 104*y[i + 1] + 114*y[i+2] - 56*y[i+3] + 11*y[i+4]
    elseif i == 2
        dy = 11*y[i-1] -20*y[i] + 6*y[i+1] + 4*y[i+2] - y[i+3]
    elseif i == Nx-1
        dy = -y[i-3] + 4*y[i-2] + 6*y[i-1] - 20*y[i] + 11*y[i+1]
    elseif i == Nx
        dy = 11*y[i-4] - 56*y[i-3] + 114*y[i-2] - 104*y[i-1] + 35*y[i]
    else
        dy = -y[i-2] + 16*y[i-1] - 30*y[i] + 16*y[i+1] - y[i+2]
    end
    dy = dy / (12 * dx^2)
    return dy
end

x = range(-10, 10, length=101)
y = @. x^2

dy = [derivative2(x, y, i) for i=1:length(x)]

@assert isapprox(dy, 2*ones(length(x)))

```

---

<div class="post-metadata">

**Author:** ![AngusRussell](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/angusrussell/32/22226_2.png) [@AngusRussell](https://discourse.julialang.org/u/AngusRussell)\
**Post date:** [October 8, 2021, 12:31pm UTC](https://discourse.julialang.org/t/best-method-for-calculating-the-second-derivative-of-1d-array/69353/7 "2021-10-08T12:31:08Z")

</div>

Thank you for your reply. The code is accurate enough for my use case, what I am now looking to do is try to optimise speed without sacrificing too much accuracy as I go on to iterate over a large range of values. This is currently taking hours to run through and the flexibility to change parameters and not wait so long for results would be great

---

<div class="post-metadata">

**Author:** ![AngusRussell](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/angusrussell/32/22226_2.png) [@AngusRussell](https://discourse.julialang.org/u/AngusRussell)\
**Post date:** [October 8, 2021, 12:32pm UTC](https://discourse.julialang.org/t/best-method-for-calculating-the-second-derivative-of-1d-array/69353/8 "2021-10-08T12:32:52Z")

</div>

Thank you very much for your reply. I will try implementing this and see how it affects the speed and accuracy of my results 🙂

---

<div class="post-metadata">

**Author:** ![Bardo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bardo/32/21601_2.png) [@Bardo](https://discourse.julialang.org/u/Bardo)\
**Post date:** [October 12, 2021, 6:29am UTC](https://discourse.julialang.org/t/best-method-for-calculating-the-second-derivative-of-1d-array/69353/9 "2021-10-12T06:29:47Z")

</div>

Further help:

- [What’s the best way to calculate a numerical derivative in MATLAB? - Stack Overflow](https://stackoverflow.com/questions/29478838/whats-the-best-way-to-calculate-a-numerical-derivative-in-matlab)
- [JuliaDiff](https://juliadiff.org/)

---

<div class="post-metadata">

**Author:** ![pitsianis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pitsianis/32/26588_2.png) [@pitsianis](https://discourse.julialang.org/u/pitsianis)\
**Post date:** [October 19, 2021, 9:34pm UTC](https://discourse.julialang.org/t/best-method-for-calculating-the-second-derivative-of-1d-array/69353/10 "2021-10-19T21:34:39Z")

</div>

To speed up the above code

- Inline the finite differences function in the loop and
- move out invariant calculations
- peel off the first two and last two loop iterations to get rid of the if statements.

---

<div class="post-metadata">

**Author:** ![AngusRussell](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/angusrussell/32/22226_2.png) [@AngusRussell](https://discourse.julialang.org/u/AngusRussell)\
**Post date:** [October 20, 2021, 4:42pm UTC](https://discourse.julialang.org/t/best-method-for-calculating-the-second-derivative-of-1d-array/69353/12 "2021-10-20T16:42:04Z")

</div>

I am quite new to Julia and programming in general so I was wondering if I have implemented your suggestions correctly. Here is my new code;

```julia-auto
# inline finite difference function
@inline function new_derivative2!(y ,i)
    # Removed if statements
    dy = -y[i-2] + 16*y[i-1] - 30*y[i] + 16*y[i+1] - y[i+2]
    dy = dy / (12 * 0.001^2)
    return dy
end

function new_QC!(E, ϵ)
    Q_C = Array{Float64, 3}(undef, 5, length(ϵ), length(ϵ)) # Define array
    dx = 0.001 # removing invariant calculation
    for k =1:5
        for j = 1:length(ϵ)
            Q_C[k,1,j] = (35*E[k,1,j] - 104*E[k,2,j] + 114*E[k,3,j] - 56*E[k,4,j] + 11*E[k,5,j]) / (12 * dx^2)
            Q_C[k,2,j] = (11*E[k,1,j] -20*E[k,2,j] + 6*E[k,3,j] + 4*E[k,4,j] - E[k,5,j]) / (12 * dx^2)
            for i=3:(length(ϵ)-2)
                Q_C[k,i,j] = new_derivative2!(E[k,:,j], i)
            end
            Q_C[k,520,j] = (-E[k,517,j] + 4*E[k,518,j] + 6*E[k,519,j] - 20*E[k,520,j] + 11*E[k,521,j]) / (12 * dx^2)
            Q_C[k,521,j] = (11*E[k,517,j] - 56*E[k,518,j] + 114*E[k,519,j] - 104*E[k,520,j] + 35*E[k,521,j]) / (12 * dx^2)
        end
    end
    return Q_C
    
end

```

Benchmarks of old and new:

## Old

```julia-auto
julia> @benchmark QC!(E, ϵ)
BenchmarkTools.Trial: 6 samples with 1 evaluation.
 Range (min … max): 925.521 ms … 1.122 s ┊ GC (min … max): 25.53% … 24.48%
 Time (median): 977.854 ms ┊ GC (median): 25.28%
 Time (mean ± σ): 995.556 ms ± 77.650 ms ┊ GC (mean ± σ): 25.16% ± 0.41%

  █ █ █ █ █ █
  █▁█▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  926 ms Histogram: frequency by time 1.12 s <

 Memory estimate: 5.51 GiB, allocs estimate: 1357207.

```

## New

```julia-auto
julia> @benchmark new_QC!(E,ϵ)
BenchmarkTools.Trial: 6 samples with 1 evaluation.
 Range (min … max): 894.819 ms … 1.077 s ┊ GC (min … max): 25.53% … 24.80%
 Time (median): 951.854 ms ┊ GC (median): 25.24%
 Time (mean ± σ): 978.050 ms ± 72.613 ms ┊ GC (mean ± σ): 25.19% ± 0.28%

  █ █ █ █ █ █
  █▁▁▁▁▁▁▁▁▁▁▁▁█▁▁█▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁█ ▁
  895 ms Histogram: frequency by time 1.08 s <

 Memory estimate: 5.47 GiB, allocs estimate: 1346787.

```

Thank you

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [October 20, 2021, 7:45pm UTC](https://discourse.julialang.org/t/best-method-for-calculating-the-second-derivative-of-1d-array/69353/13 "2021-10-20T19:45:06Z")

</div>

One last thing is that `dy = dy / (12 * 0.001^2)` will be a good bit slower than `dy = dy * (1/ (12 * 0.001^2))` since the second version will use a multiplication by a constant and the first will use a division.