I’m learning how to solve PDEs numerically.

The code I wrote in Julia 1.2.0 is as follows:

```
using PyCall
so = pyimport("scipy.optimize")
using PyPlot
I = im
h_s = 0.1
h_t = 0.02
T = 100.0
L_x = 30.0
L_y = 30.0
ddx(u) = (circshift(u, (1,0)) - circshift(u, (-1,0)))/(2*h_s)
ddy(u) = (circshift(u, (0,1)) - circshift(u, (0,-1)))/(2*h_s)
curl(v_x, v_y) = ddx(v_y) - ddy(v_x)
div(v_x, v_y) = ddx(v_x) + ddy(v_y)
laplacian(u) = (circshift(u, (1,0)) + circshift(u, (-1,0)) + circshift(u, (0, 1)) + circshift(u, (0,-1)) - 4*u)/(h_s^2)
# Define the mesh
mesh_x = 0:h_s:L_x
mesh_y = 0:h_s:L_y
mesh_T = 0:h_t:T
# Indices for start and end points of the domain
x_i = 1
x_f = length(mesh_x)
y_i = 1
y_f = length(mesh_x)
# Declare arrays storing the result
gridSize = (length(mesh_x), length(mesh_y))
u = [zeros(gridSize) for t=mesh_T]
# Initial condition
u[1] = [exp(-((x-2)^2+(y-2)^2)) for x=mesh_x, y=mesh_y]
function BC(U)
U[x_i,:] .= 0.0
U[x_f,:] .= 0.0
U[:,y_i] .= 0.0
U[:,y_f] .= 0.0
return U
end
function F(U)
U = reshape(U, (x_f, y_f))
return vec(BC(laplacian(U)))
end
function Crank_Nicolson(U)
U = vec(BC(U))
function f(u_next)
return u_next - U - h_t/2*(F(U)+F(u_next))
end
sol = so.root(f, U, method="krylov")
return reshape(sol["x"], (x_f, y_f))
end
for t in 1:length(mesh_T)-1
u[t+1] = Crank_Nicolson(u[t])
end
# imshow(u[end])
# show()
```

I have a python version of this code.

The time of running for julia is about 6m30s, and 4m40 for python code.

I noticed there is a “warmup” period for julia and it’s about 20s.

Even considering this, it seems python is still faster than julia.

The most time-consuming part of this program is solving the time evolution.

Crank-Nicolson is an implicit method so there is an nonlinear equation to be solved.

For this part, I use the scipy.optimize.root function for both julia and python version.

Thus I would expect both versions should spend similar time,

but the result is julia spent more than 25% time than python.

Does anyone have ideas on this?

Thank you.