2015-12-08 9:10 GMT-05:00 博陈 <[email protected]>:
> 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
`const` in local scope is basically no-op
> 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
Don't default functions in local scope.
> 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)
Also apply to code below, if you really care about performance (e.g.
your loop below), pre-allocate the buffer and use in place operations
instead.
>
> #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])
Very minor, there's `abs2`.
> #show()
I usually also prefer to put these plotting functions outside the
function I benchmark so that they don't mess up timing or @code_* and
I can check the result at the same time.
> end
>
> @time ground()
>
