# Is there a central difference/gradient function somewhere?

**URL:** <https://discourse.julialang.org/t/is-there-a-central-difference-gradient-function-somewhere/19454>\
**Category:** General Usage\
**Created:** [January 9, 2019, 10:30pm UTC](https://discourse.julialang.org/t/is-there-a-central-difference-gradient-function-somewhere/19454 "2019-01-09T22:30:27Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![weech](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weech/32/6717_2.png) [@weech](https://discourse.julialang.org/u/weech)\
**Post date:** [January 9, 2019, 10:30pm UTC](https://discourse.julialang.org/t/is-there-a-central-difference-gradient-function-somewhere/19454/1 "2019-01-09T22:30:27Z")

</div>

When I write in Python, I often use the gradient function from NumPy, and, when I write in NCL, I often use the center\_finite\_diff\_n function. However, it looks like Julia does not have a function that has a behavior similar to those. I tried searching for an equivalent function among the 3rd-party packages, but couldn’t find anything that takes arrays, only functions. Does anyone know of any package with a function that behaves like NumPy’s gradient, NCL’s center\_finite\_diff\_n, or Matlab’s gradient?

---

<div class="post-metadata">

**Author:** ![Mattriks](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mattriks/32/351_2.png) [@Mattriks](https://discourse.julialang.org/u/Mattriks)\
**Post date:** [January 10, 2019, 8:55am UTC](https://discourse.julialang.org/t/is-there-a-central-difference-gradient-function-somewhere/19454/2 "2019-01-10T08:55:30Z")

</div>

There is a function `gradvecfield()` in [`CoupledFields.jl`](https://github.com/Mattriks/CoupledFields.jl). It’s used by [`Gadfly.jl`](https://github.com/GiovineItalia/Gadfly.jl/) to calculate the gradient vector field for [`Geom.vectorfield`](http://gadflyjl.org/stable/gallery/geometries/#%5BGeom.vectorfield%5D(@ref)-1). `gradvecfield()` works like this:

```julia
using CoupledFields
kernelpars = GaussianKP(X)
∇g = gradvecfield([a b], X, Y, kernelpars)

```

where `a` is a smoothness parameter, `b` is a ridge parameter, `X` and `Y` are matrices, and `Y = g(X)`.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [January 10, 2019, 8:59am UTC](https://discourse.julialang.org/t/is-there-a-central-difference-gradient-function-somewhere/19454/3 "2019-01-10T08:59:11Z")

</div>

something like this?

```julia
function centraldiff(v::AbstractMatrix)
    dv = diff(v)/2
    a = [dv[[1],:];dv]
    a .+= [dv;dv[[end],:]]
    a
end

function centraldiff(v::AbstractVector)
    dv = diff(v)/2
    a = [dv[1];dv]
    a .+= [dv;dv[end]]
    a
end

```

These methods are of course not optimized for performance, writing the loop manually would likely improve performance and reduce memory allocations.

Note, these functions do not reduce the length of the input array (as `diff` does), by copying the first and last `diff` elements.

---

<div class="post-metadata">

**Author:** ![mauro3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mauro3/32/292_2.png) [@mauro3](https://discourse.julialang.org/u/mauro3)\
**Post date:** [January 10, 2019, 12:18pm UTC](https://discourse.julialang.org/t/is-there-a-central-difference-gradient-function-somewhere/19454/4 "2019-01-10T12:18:38Z")

</div>

There is [https://github.com/JuliaDiffEq/DiffEqOperators.j](https://github.com/JuliaDiffEq/DiffEqOperators.j); intro docs [https://github.com/JuliaDiffEq/DiffEqOperators.jl/blob/master/docs/DiffEqOperators.md](https://github.com/JuliaDiffEq/DiffEqOperators.jl/blob/master/docs/DiffEqOperators.md).

---

<div class="post-metadata">

**Author:** ![weech](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weech/32/6717_2.png) [@weech](https://discourse.julialang.org/u/weech)\
**Post date:** [January 10, 2019, 5:33pm UTC](https://discourse.julialang.org/t/is-there-a-central-difference-gradient-function-somewhere/19454/5 "2019-01-10T17:33:27Z")

</div>

While I appreciate your responses, all the functions you pointed out only work on 2d arrays. I put together a function that works on 4d arrays and operates on the 3rd dimension with irregular spacing on the derivative dimension (aka the use case I needed yesterday).

```julia
function partialp(arr, coord)

    nx, ny, np, nt = size(arr)
    out = similar(arr)
    dcoord = diff(coord)

    # Forward difference at bottom
    dp = dcoord[1]
    for t in 1:nt, y in 1:ny, x in 1:nx
        out[x, y, 1, t] = (arr[x, y, 2, t] - arr[x, y, 1, t]) / dp
    end

    # Central difference in interior using numpy method
    for p in 2:np-1
        dp1 = dcoord[p-1]
        dp2 = dcoord[p]
        a = -dp2 / (dp1 * (dp1 + dp2))
        b = (dp2 - dp1) / (dp1 * dp2)
        c = dp1 / (dp2 * (dp1 + dp2))
        for t in 1:nt, y in 1:ny, x in 1:nx
            out[x, y, p, t] = a*arr[x, y, p-1, t] + b*arr[x, y, p, t] + c*arr[x, y, p+1, t]
        end
    end

    # Backwards difference at top
    dp = dcoord[end]
    for t in 1:nt, y in 1:ny, x in 1:nx
        out[x, y, end, t] = (arr[x, y, end, t] - arr[x, y, end-1, t]) / dp
    end

    return out
end

```

It produces the same results as numpy.gradient, is faster than using PyCall, and doesn’t make many extraneous allocations. However, it doesn’t have any of the checks that a function would need to be used generally, and it’s just as limited as the matrix examples you gave, but it would be a good jumping off point for anyone trying to implement a central difference calculation in Julia.

Edit: Swapped t and x in loops to take advantage of Julia being column major.

---

<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:** [January 10, 2019, 5:59pm UTC](https://discourse.julialang.org/t/is-there-a-central-difference-gradient-function-somewhere/19454/6 "2019-01-10T17:59:23Z")

</div>

> [@weech](#):
>
> Does anyone know of any package with a function that behaves like NumPy’s gradient, NCL’s center\_finite\_diff\_n, or Matlab’s gradient?

Used to be in Base, but was removed (see [`gradient()`: remove from Base, then deprecate · Issue #16113 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/issues/16113)). That implementation only ever supported 1d arrays, though.

Something could be added to [GitHub - MatlabCompat/MatlabCompat.jl: Source of MatlabCompat.jl](https://github.com/MatlabCompat/MatlabCompat.jl)

---

<div class="post-metadata">

**Author:** ![weech](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weech/32/6717_2.png) [@weech](https://discourse.julialang.org/u/weech)\
**Post date:** [January 10, 2019, 6:20pm UTC](https://discourse.julialang.org/t/is-there-a-central-difference-gradient-function-somewhere/19454/7 "2019-01-10T18:20:35Z")

</div>

Any general implementation of the central difference would be best based off of NumPy’s, since Matlab’s gradient doesn’t handle non-uniform spacing. Reading through NumPy’s implementation, it looks like it’s relatively simple to port to Julia. That might be a good weekend project.

---

<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:** [January 10, 2019, 6:22pm UTC](https://discourse.julialang.org/t/is-there-a-central-difference-gradient-function-somewhere/19454/8 "2019-01-10T18:22:50Z")

</div>

> [@weech](#):
>
> When I write in Python, I often use the `gradient` function from NumPy

By the way, I’m curious: what do you use this function for? Plotting?

---

<div class="post-metadata">

**Author:** ![Mattriks](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mattriks/32/351_2.png) [@Mattriks](https://discourse.julialang.org/u/Mattriks)\
**Post date:** [January 10, 2019, 6:27pm UTC](https://discourse.julialang.org/t/is-there-a-central-difference-gradient-function-somewhere/19454/9 "2019-01-10T18:27:12Z")

</div>

@weech To clarify `CoupledFields.gradvecfield` works on arrays of any dimension, and for irregular spaced points (its implementation in Gadfly is only for 2d arrays). Say X is a n\times p matrix and Y is a n\times q matrix, then `CoupledFields.gradvecfield([a b], X, Y, kernelpars)` returns n gradient matrices (of size p\times q), for the n function points in `X`.

---

<div class="post-metadata">

**Author:** ![weech](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weech/32/6717_2.png) [@weech](https://discourse.julialang.org/u/weech)\
**Post date:** [January 10, 2019, 6:30pm UTC](https://discourse.julialang.org/t/is-there-a-central-difference-gradient-function-somewhere/19454/10 "2019-01-10T18:30:32Z")

</div>

> [@stevengj](#):
>
> By the way, I’m curious: what do you use this function for? Plotting?

Nope, in today’s case I need to calculate isobaric [potential vorticity](https://en.wikipedia.org/wiki/Potential_vorticity). I’m calculating the zonal and meridional derivatives using Spherepack, but I needed the central difference for the vertical derivative.

---

<div class="post-metadata">

**Author:** ![weech](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weech/32/6717_2.png) [@weech](https://discourse.julialang.org/u/weech)\
**Post date:** [January 10, 2019, 6:36pm UTC](https://discourse.julialang.org/t/is-there-a-central-difference-gradient-function-somewhere/19454/11 "2019-01-10T18:36:35Z")

</div>

> [@Mattriks](#):
>
> @weech To clarify `CoupledFields.gradvecfield` works on arrays of any dimension, and for irregular spaced points (its implementation in Gadfly is only for 2d arrays). Say X is a n\times p matrix and Y is a n\times q matrix, then `CoupledFields.gradvecfield([a b], X, Y, kernelpars)` returns n gradient matrices (of size p\times q), for the n function points in `X`.

Oh, OK, I was confused because the signature of that function was

```julia
function gradvecfield(par::Array{Float64}, X::T, Y::T, kpars::KernelParameters ) where T<:Matrix{Float64}

```

and got confused by the Matrix typing. It’s still several levels beyond what I know about mathematics, so it’d take me a long time to figure out how to use it.

---

<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:** [January 10, 2019, 6:37pm UTC](https://discourse.julialang.org/t/is-there-a-central-difference-gradient-function-somewhere/19454/12 "2019-01-10T18:37:25Z")

</div>

> [@weech](#):
>
> Nope, in today’s case I need to calculate isobaric [potential vorticity](https://en.wikipedia.org/wiki/Potential_vorticity).

As is often the case when translating vectorized code from other languages to Julia, there is probably a better way to do it because Julia doesn’t force you to vectorize your algorithms for performance.

A calculation of potential vorticity based on `gradient` calls [makes a lot of passes over the array(s)](https://julialang.org/blog/2017/01/moredots#why-vectorized-code-is-not-as-fast-as-it-could-be) and allocates several temporary arrays, all to compute a scalar result at the end. Instead, in Julia, you could do it all in a single (nested) loop over the data (one pass, no temporaries). I wouldn’t be surprised if a properly coded loop were an order of magnitude faster than a method based on `gradient` calls.

---

<div class="post-metadata">

**Author:** ![weech](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weech/32/6717_2.png) [@weech](https://discourse.julialang.org/u/weech)\
**Post date:** [January 10, 2019, 6:51pm UTC](https://discourse.julialang.org/t/is-there-a-central-difference-gradient-function-somewhere/19454/13 "2019-01-10T18:51:47Z")

</div>

> [@stevengj](#):
>
> I wouldn’t be surprised if a properly coded loop were an order of magnitude faster than a method based on `gradient` calls.

Hmm, that might not be a bad idea. I’ve been using a certain formula (Bluestein’s Synoptic-Dynamic Meteorology in Midlatitudes. Eq 4.5.93) that uses x-y-p coordinates, and it might be worth deriving something in ϕ-θ-p coordinates so I can use it with my data without the spherical transforms.