Solving-Parabolic-PDEs-with-finite-differences-methods

Solution of the heat/diffusion and time-dependent Schrodinger equations with Crank-Nicolson and Euler methods

By: Alexandros Stratoudakis

Date: Spring 2024

---

The project aims to solve parabolic PDEs, namely the diffusion and Schrodinger equations in 1d and 2d.

This file is instructions to the basic usage of the code. Most classes and functions used, have a thorough doc-sheet as a comment where the parameters and returns are thoroughly described. If someone is interested in such details, he hasn't but to open the files and read the comments on the definitions of functions and classes.

---

Code Structure:

1. The bulk of the code is in diffusion.py and schrodinger.py. There, simulations are defined as python classes and solvers are called as methods of these classes.

2. There are four simulation classes: Diffusion_1d, Diffusion_2d, Schrodinger_1d and Schrodinger_2d. All, except 2D diffusion, use the Crank-Nicolson method to solve the equations using the method .Crank_Nicolson. For 2D diffusion, Euler's method is used instead with .Euler .

3. The utilities.py script, contains useful function such as potentials, initial distributions and wave-packets.

4. All scripts that start with "run" (e.g. run_schrodinger_2D.py) are non-essential. They are examples of runfiles that use the code for some cases. With minor modifications to these, one can explore the full capabilities of the code.

---

How to run the code:

As mentioned before the example runfiles provided can be a good starting point for any application. However, I will outline the process in case someone wants to make his own.

1. It is strongly advised to open the parent (to the files) directory "as a project" with either Spyder (suggested) or PyCharm.

2. Import the class of the desired solver along with anything you need from the utilities script.

3. Instantiate the simulation, e.g. sim = Schrodinger_1d( L = 50., x0=0, T = 10, Nx = 2501, Nt = 10001 )

4. Use the solver method, e.g. u = sim.Crank_Nicolson(lambda x: gaussian_wave_packet_1d(x, lamda=8, center = 10, k0 = 10) , Left_BC = 0., Right_BC = 0., BC_type = 'dirichlet', normalize_input = True, V=0)

5. Plot and/or save results e.g. with matplotlib.

The above example will simulate the propagation of a Gaussian wave-packet.

Note that in many cases, especially in 2D, the program can take significant time to run.

---

Examples

1. 1D alpha radiation

gaussian_wavepacket_1d_alpha_E.78.125Vmin.-50Vmax.91.09.mp4

2. 1D scattering by a step-potential

gaussian_wavepacket_1d_scattering_heaviside_E.12.5V0.15.0.mp4

3. Particle in a box time-evolution

particle_in_a_box_ground_state.mp4

particle_in_a_box_psi1+psi2.mp4

4. 2D diffusion

diff2D_donut_dirichlet.mp4

5. 2D double-slit experiment

double_slit_E50.0_w0.2_th0.2_d0.8_dx0.05_dt0.0005.mp4

References

1. Trachanas, S. (2018). An Introduction to Quantum Physics: A First Course for Physicists, Chemists, Materials Scientists, and Engineers. Wiley-VCH.
2. Crank, J., & Nicolson, P. (1947). A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type. Mathematical Proceedings of the Cambridge Philosophical Society, 43(1), 50-67. doi:10.1017/S0305004100023197
3. Euler, L. (1768). Institutiones calculi differentialis cum eius usu in analysi finitorum ac doctrina serierum. IMPENSIS ACADEMIAE IMPERIALIS SCIENTIARUM. Retrieved from https://archive.org/details/institutionescal00eule
4. Excellent repo by artmenlope (https://github.com/artmenlope/double-slit-2d-schrodinger) that partially inspired my 2D Schrodinger Crank-Nicolson scheme.