Are you sure this version works? I get Δϕ == 1.0 in all iterations, so I
can't see how it would ever exit the while loop...
On Monday, December 14, 2015 at 4:35:43 AM UTC+1, 博陈 wrote:
>
> 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
>
>
