This is my first julia code, I am happy it did the right thing, but compared with the Matlab code that did the same thing, it runs so slowly. The Matlab code takes about 90s, but the julia below takes 130s.
Potential(x::Float64, y::Float64) = -2./sqrt(x^2+1.) -2./sqrt(y^2+1.) + 1./sqrt((x-y)^2+1.) Ko(x::Float64,y::Float64) = x^2+y^2 function ground() R = 320. Ks = 2^11 x = linspace(-R, R, Ks+1) dx = 2R/Ks dt = 0.1 pm = 2π/2dx px = circshift( linspace(-pm, pm, Ks+1), round(Int64, Ks/2) ) FFTW.set_num_threads(8) V = Float64[Potential(ix, iy) for ix in x, iy in x] ko = Float64[e^(-dt/4* (Ko(ipx,ipy))) for ipx in px, ipy in px] ko2 = ko.^2 vo = Float64[e^(-dt* (Potential(ix,iy))) for ix in x, iy in x] ϕ = Array(Complex128,(Ks+1,Ks+1)) Pl = plan_fft!(ϕ; flags=FFTW.MEASURE) ϕ = V.*complex(1.) normphi = sqrt(sumabs2(ϕ))*dx ϕ /= normphi ϕ = Pl*ϕ ϕ .*= ko ϕ = Pl\ϕ Δϕ = 1. nstep = 0 print("start loop") while Δϕ > 1.e-15 ϕ′ = ϕ ϕ .*= vo normphi = sqrt(sumabs2(ϕ))*dx ϕ /= normphi ϕ = Pl*ϕ ϕ .*= ko2 ϕ = Pl\ϕ # if check Δϕ for every step, 35s is needed. if nstep>500 Δϕ = maxabs(ϕ-ϕ′) end nstep += 1 if mod(nstep,200)==0 print(nstep," ") print("Δϕ=",Δϕ," ") end end ϕ .*= vo ϕ = Pl*ϕ ϕ .*= ko ϕ = Pl\ϕ normphi = sqrt(sumabs2(ϕ))*dx ϕ /= normphi end