search

maroba / findiff

Python package for numerical derivatives and partial differential equations in any number of dimensions.
MIT License
438 stars 61 forks source link

Using findiff to solve the stationary Schrödinger equation #49

Closed yyanghy closed 2 years ago

yyanghy commented 2 years ago

When trying to solve the stationary Schrödinger equation following the example at https://findiff.readthedocs.io/en/latest/source/examples-matrix.html, the code doesn't work.

x1_list=np.linspace(-8 pi,8 pi,201) x2_list=np.linspace(-0.5 pi,1.5 pi,201) dx=x1_list[1]-x1_list[0] dy=x2_list[1]-x2_list[0] x,y=np.meshgrid(x1_list,x2_list)

u=(-2 20 np.cos(y) np.cos(x-pi) + 0.03 x**2)

Hm=(-15 FinDiff(0,dx,2)-0.04 FinDiff(1,dy,2)).matrix(u.shape)+u.reshape(-1) energies, states = eigs( Hm, k=5, which='SR' ) bit0=states[:,0].reshape(201,201)

There is a memory error: MemoryError: Unable to allocate 12.2 GiB for an array with shape (40401, 40401) and data type float64

maroba commented 2 years ago

The matrix representation of a FinDiff operator is a sparse matrix, and so is your kinetic term in the Hamiltonian. However, when you add the potential by saying +u.reshape(-1), you are turning the whole Hamiltonian into a dense matrix representation (although most values are zero). This eats up all your memory.

What you need to do: declare your potential also as a sparse matrix. Someone has written an article at Medium.com for this use case.

yyanghy commented 2 years ago

Thank you so much for your reply!