I am pretty new to Julia and ultimately I want to transfer my python code for solving a nonlinear PDE into Julia (and see whether I can get the speedup I hope).
The PDE I want to solve is similar to a nonlinear heat equation. In python I have a function which computes dh/dt. I use the method of lines to discretize the vector h(x,t) in spatial domain, and then use a stiff ODE solver scipy.integrate.solve_ivp(… method=‘BDF’) to solve the PDE. I am trying to adopt a similar approach in Julia with the DifferentialEquations.jl package.
The first step I need to do is to figure out a way to implement numerical differentiation using finite difference (I need first and third derivative, with either 2nd or 4th order accuracy) with periodic boundary condition applied to both ends of the vector. I have the following four versions of codes:
function ndiff1(u::AbstractVector, n::Integer, order::Integer, dx::Float64)
Nx = length(u)
if n==1
if order==2
du = 0.5*(circshift(u,1)  circshift(u,1))
elseif order==4
du = (circshift(u,2)  circshift(u,2) + 8*(circshift(u,1)  circshift(u,1)))/12.0
end
elseif n==3
# similar to above with just different stencils
end
dudx = du/dx^n
return dudx
end
function ndiff2(u::AbstractVector, n::Integer, order::Integer, dx::Float64)
Nx = length(u)
dudx = similar(u)
if n==1
if order==2
for i = 2:Nx1
dudx[i] = 0.5*(u[i+1]  u[i1])/dx
end
dudx[1] = 0.5*(u[2]u[end])/dx
dudx[end] = 0.5*(u[1]u[end1])/dx
elseif order==4
# similar to above with different stencils
end
elseif n==3
# similar to above with different stencils
end
return dudx
end
using DiffEqOperators
const periodic_bc = PeriodicBC(Float64)
const D1_O2 = CenteredDifference(1, 2, dx, Nx)*periodic_bc
const D1_O4 = CenteredDifference(1, 4, dx, Nx)*periodic_bc
const D3_O2 = CenteredDifference(3, 2, dx, Nx)*periodic_bc
const D3_O4 = CenteredDifference(3, 2, dx, Nx)*periodic_bc
function ndiff3(u::AbstractVector, n::Integer, order::Integer)
if n==1
if order==2
dudx = D1_O2*u
elseif order==4
dudx = D1_O4*u
end
elseif n==3
if order==2
dudx = D3_O2*u
elseif order==4
dudx = D3_O4*u
end
end
return dudx
end
function ndiff4(u::AbstractVector, n::Integer, order::Integer, dx::Float64)
Nx = length(u)
dudx = CenteredDifference(n, order, dx, Nx)*periodic_bc*u
return dudx
end
When I used @btime to benchmark the above four functions on an array of size 256 (typical for my problem), ndiff3
is a bit faster than ndiff2
, which is a bit faster than ndiff1
, and ndiff4
is close to an order of magnitude slower than the rest. My questions are:

Are there anything noticeable areas for improvement for
ndiff1
andndiff2
? That’s similar to how I implemented the finite difference in python (except I try to ‘vectorize’ it in python using numpy without the for loop forndiff2
) and that’s already close to 2 orders of magnitude of speedup which I’m pretty happy about, but I’m not sure if I there is anything significant I’m missing in Julia. 
I try to use the
DiffEqOperators
package forndiff3
andndiff4
so I don’t have to write everything explicitly from scratch myself. However I’m not sure what’s the best way to implementCenteredDifference(n, order, dx, Nx)*periodic_bc
. If I place them inside the function as inndiff4
, does Julia need to initialize it every time the function is called (and hence leads to slower speed)? If I place them outside the function as I did forndiff3
, then obviously it’s very inflexible as I will need to define the problem sizeNx
beforehand.
any comment/help is appreciated. thanks