I get it, thanks for pointing that out.
On Monday, December 14, 2015 at 11:38:14 AM UTC+1, FQ wrote:
>
> The code does exit at some point, returning a 2049x2049 complex matrix.
> note that delta phi is reassigned inside the loop at some point.
>
> Am 14.12.2015 um 10:31 schrieb Patrick Kofod Mogensen:
> > 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
> >
>
>
