The .* and ./ operators are slow. This runs in about 8 seconds. I get a bunch of Nans as a result so maybe I messed up somewhere but you can at least see the same principles:
https://gist.github.com/KristofferC/a7d61aa7cb9ea7e18d3e Note that I used the Devectorize package. On Tuesday, December 8, 2015 at 4:17:41 PM UTC+1, 博陈 wrote: > > function ground() > #using PyPlot > > print("start") > const R = 320. > const Ks = 2^11 > const x = linspace(-R, R, Ks+1) > const dx = 2R/Ks > const y = x > const dt = 0.1 > const pm = 2π/2dx > px = circshift( linspace(-pm, pm, Ks+1), round(Int64, Ks/2) ) > py = px > FFTW.set_num_threads(8) > #Potential(x::Float64, y::Float64) = -2/sqrt((x-2.)^2+1) > -2/sqrt((y-2)^2+1) -1/sqrt((x+2.)^2+1) -1/sqrt((y+2)^2+1) + > 1/sqrt((x-y)^2+1) > 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 > V = Float64[Potential(ix, iy) for ix in x, iy in y] > ϕ = V.*(1.+0im) > ϕ /= maximum(abs(ϕ)) > ko = Float64[e^(-dt/4* (Ko(ipx,ipy))) for ipx in px, ipy in py] > ko2 = ko.^2 > #ϕ = ground(x,V) > > #function ground(x, V) > > ϕ = fft(ϕ) > ϕ .*= ko > ϕ = ifft(ϕ) > Δϕ = 1. > nstep = 0 > while Δϕ > 1.e-15 > ϕ′ = ϕ > ϕ .*= e.^(-dt*V) > ϕ /= maximum(abs(ϕ)) > ϕ = fft(ϕ) > ϕ .*= ko2 > ϕ = ifft(ϕ) > > Δϕ = maximum(abs(ϕ-ϕ′)) > nstep += 1 > if mod(nstep,200)==0 > print(nstep) > print("Δϕ=",Δϕ) > end > end > ϕ .*= exp(-dt*V) > ϕ = fft(ϕ) > ϕ .*= ko > ϕ = ifft(ϕ) > ϕ /= maximum(abs(ϕ)) > > #imshow(abs(ϕ)^2, interpolation="none", extent=[-R, R ,-R ,R]) > #show() > end > > @time ground() > >
