Hi, Steven. This is how I'm building the matrix:
#======================#
# OPERATING THE MATRIX #
#======================#
# --Creating HAMILTONIAN container.
M = spzeros(n*n, n*n)#--Creating the SPARSE matrix (n is
the number of points per row).
# --Working kinetic energy into matrix M (using 2D Laplacian scheme)...
# --Central part:
trmC = -4#--el termino central de la matriz.
M[1,1] = trmC
for i = 2:n*n M[i,i] = trmC; M[i-1,i] = 1 end#--diags. central and 1
above.
for i = 1:(n*n)-1 M[i+1,i] = 1 end#--diag. below.
#--Exterior diagonals (2 in 1 loop).
for i = 1:n*(n-1) M[i,i+n] = 1; M[i+n,i] = 1 end
#--"poking" zeroes into the matrix to make "holes."
for i = 1:n*(n-1) if i%n==0 M[i+1,i]=0;M[i,i+1] = 0 end end
After that I multiply M by the appropriate constants to get the full
kinetic energy. Then I use the routine I posted previously so as to add the
potential, all in the same matrix.
I assume my matrix creation scheme may not be the most efficient, but it
did get the job done (up til now).
Cheers.
El miércoles, 2 de noviembre de 2016, 3:43:48 (UTC-7), Alejandro
Castellanos escribió:
>
> Hello.
>
> I am working with a pentadiagonal sparse matrix that represents a 2D
> Schrodinger's time-independent equation. I first work the laplacian
> expressed in Finite Differences form and then I apply the potential on the
> same matrix.
>
> So far I've been able to validate my results for both an electron in a box
> as well as a harmonic oscillator, but when I change to the following
> potential of a dipole, Julia pretty much quits on me when I try to obtain
> the eigenvalues and eigenvectors:
>
> O = [round(L/2)-hx round(L/2)-hy]# --el ORIGEN (centro -- x,y)
> del potencial.
> Eps_o = 8.854187817e10-12# --F*m^-1
> C = 1/(4*pi*Eps_o)
> D = 1e-21#C*m^2/s# --Debyes)
> pe = 1.8*D
> *P(X,Y) = -(C)*pe*(Y/X)*(1/( (X)^2 + (Y)^2 ) )*# --How the
> potential gets described.
>
> #--I'm aware there's singularities in the potential.
> #--and here's how I apply the potential to my sparse matrix.
>
> Vi = Float64[]# --container for the potential.
> for j=Y for i=X push!(Vi,P(i,j)) end end@ --applying the potential.
>
>
> I use this command: *l, v = eigs(M,nev=15,which = :SM ,ritzvec=true)*
>
> My problem seems to be that there's an error that I can't get past:
>
> ERROR: LoadError: ArgumentError: matrix has one or more zero pivots
>>>
>>> in #ldltfact!#10(::Float64, ::Function,
>>>> ::Base.SparseArrays.CHOLMOD.Factor{Float64},
>>>> ::Base.SparseArrays.CHOLMOD.Sparse{Float64}) at ./sparse/cholmod.jl:1350
>>>
>>> in (::Base.LinAlg.#kw##ldltfact!)(::Array{Any,1},
>>>> ::Base.LinAlg.#ldltfact!, ::Base.SparseArrays.CHOLMOD.Factor{Float64},
>>>> ::Base.SparseArrays.CHOLMOD.Sparse{Float64}) at ./<missing>:0
>>>
>>> in #ldltfact#12(::Float64, ::Array{Int64,1}, ::Function,
>>>> ::Base.SparseArrays.CHOLMOD.Sparse{Float64}) at ./sparse/cholmod.jl:1386
>>>
>>> in #ldltfact#13(::Array{Any,1}, ::Function,
>>>> ::Hermitian{Float64,SparseMatrixCSC{Float64,Int64}}) at
>>>> ./sparse/cholmod.jl:1426
>>>
>>> in factorize(::SparseMatrixCSC{Float64,Int64}) at ./sparse/linalg.jl:897
>>>
>>> in #_eigs#62(::Int64, ::Int64, ::Symbol, ::Float64, ::Int64, ::Void,
>>>> ::Array{Float64,1}, ::Bool, ::Base.LinAlg.#_eigs,
>>>> ::SparseMatrixCSC{Float64,Int64}, ::UniformScaling{Int64}) at
>>>> ./linalg/arnoldi.jl:251
>>>
>>> in (::Base.LinAlg.#kw##_eigs)(::Array{Any,1}, ::Base.LinAlg.#_eigs,
>>>> ::SparseMatrixCSC{Float64,Int64}, ::UniformScaling{Int64}) at ./<missing>:0
>>>
>>> in #eigs#55(::Array{Any,1}, ::Function,
>>>> ::SparseMatrixCSC{Float64,Int64}, ::UniformScaling{Int64}) at
>>>> ./linalg/arnoldi.jl:78
>>>
>>> in (::Base.LinAlg.#kw##eigs)(::Array{Any,1}, ::Base.LinAlg.#eigs,
>>>> ::SparseMatrixCSC{Float64,Int64}, ::UniformScaling{Int64}) at ./<missing>:0
>>>
>>> in #eigs#54(::Array{Any,1}, ::Function,
>>>> ::SparseMatrixCSC{Float64,Int64}) at ./linalg/arnoldi.jl:77
>>>
>>> in (::Base.LinAlg.#kw##eigs)(::Array{Any,1}, ::Base.LinAlg.#eigs,
>>>> ::SparseMatrixCSC{Float64,Int64}) at ./<missing>:0
>>>
>>> in include_from_node1(::String) at ./loading.jl:488
>>>
>>> while loading /home/alejandro/Desktop/ACAD/PROG/ACADEMIC_PROGRAMMING/FDM
>>>> (Finite_Difference_Method)/2D/SCHROD_DIP_2D/SCRATCH-2D-SI-DIPOLO2.jl, in
>>>> expression starting on line 106
>>>
>>>
>>> My question is, is there a way to work around it, or, am I completely
> screwed?
>
> Thanks so much in advance.
>
