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 = [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?