# Numerical Solution of Schrödinger Equation by Crank–Nicolson Method

Authors: Amin Khan; Muhammad Ahsan; Ebenezer Bonyah; Rashid Jan; Muhammad Nisar; Abdel-Haleem Abdel-Aty; Ibrahim S. Yahia
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In this study, we implemented the well-known Crank–Nicolson scheme for the numerical solution of Schrödinger equation. The numerical results converge to the exact solution because the Crank–Nicolson scheme is unconditionally stable and accurate. We have compared the results for different parameters with analytical solution, and it is found that the Crank–Nicolson scheme is suitable for the numerical solution of Schrödinger equations. Three different problems are included to verify the accuracy, stability, and capability of the Crank–Nicolson scheme.


---

## Body

## 1. Introduction

Differential equations are widely used in many fields and are useful tools for analyzing dynamic phenomena. They also play a very important role in description of physical system into mathematical form. Diversities of physical system such as the flow of current through a conductor, the growth rate of the population of some vintage (species), the motion of a missiles, potential of uniform electric and magnetic field, the operation of the mixture, the spread of different diseases, and the pull of gravitational field of Earth are modeled by differential equation.The model equations which are depending only on one independent variable are helpful in science and engineering field particularly in the field of classical mechanics. Some of the physical problems that can be modeled by using ordinary differential equations are Newton’s law of cooling, the amount of glucose absorbs by the body, spread of epidemics, population growth, simple harmonic motion, Newton’s second law of motion, and the motion of free falling bodies [1–3].The physical system depending on more than one independent variable is commonly used to model the complex dynamic system in applied science such as biology and finance. Many experts frequently used partial differential equations to model various dynamic systems of their interest based on their knowledge and understanding. The parameters in the equation often have interesting scientific implications, but their values are frequently unknown and must be estimated from dynamic system having measurement errors [4, 5]. These types of equation have no analytical solution and can only be solved with numerical methods.The Schrödinger equation is one of the basic equation of the quantum theory. It is also the fundamental equation of quantum mechanics, which deals with submicroscopic events [6]. Erwin Schrödinger, an Austrian physicist, developed this equation in 1926. The Schrödinger equation is an essential to quantum mechanics as Newton’s laws of motion to classical mechanics large-scale phenomena. The study of this equation plays an important role in modern physics [7]. It also describes how external forces affect the structure of probability waves (see [8]) that govern the motion of tiny particles. Schrödinger demonstrated the equation’s correctness by applying it to the hydrogen atom, which accurately predicted many of its properties. It has also applications in atomic, nuclear, and solid-state physics.The Schrödinger equation has two forms, one of which explicitly includes time and so describes how a particle’s wave function would evolve over time. The wave function is often referred to as the time-dependent Schrödinger wave equation since it acts like a wave in general. The other is the steady-state Schrödinger equation. The time-dependent equation is used to describe progressive waves, which are applicable to free particle motion, whereas Schrödinger’s time-independent wave equation is used to describe standing waves. The particle’s potential energy does not always depend on time, and the potential energy is only a function of position. In such instances, the particle’s behavior is described by Schrödinger’s time-independent wave equation. [9, 10]. In general, the solutions to the time-dependent Schrödinger equation will explain the particle’s dynamical behavior in a manner comparable to Newton’s equation F=ma in classical physics [11, 12]. The nonlinear Schrödinger equation also appears in applied physics having important applications in various fields including plasma physics, nonlinear optics, and water waves [13–15]. The nonlinear Schrödinger equation has been investigated by Yajima and Outi [16], Karpman and Krushkal [17], Satsuma and Yajima [18], Hardin and Tappert [19, 20], and Tappert [21] using various techniques. Mathematically, the Schrödinger equation can be written as(1)ιℏ∂ux,t∂t=−ℏ22m∇2ux,t+Vx,tux,t,where ι=−1 represents the imaginary part, ℏ is Planck’s constant having unit in joule second, m is the mass of the particle, ux,t is the wave function define over space and time, ∇2 is the Laplacian operator, and Vx,t is the potential energy influencing the particle. This equation is the important equation of the quantum mechanics related to matter in terms of the wave-like properties of particles in a field.The numerical treatment is an alternative approach to get the required solution. Some of the recent numerical contribution related to Schrödinger types equation are studied in [6, 22] by a stable Haar wavelet collocation method (HWCM). In [23], the cubic Schro..dinger equation has been converted into system of PDEs by using ux,t=rx,t+isx,t that gives separate equations for the real and imaginary part which are then solved by HWCM. The two-dimensional Schro..dinger equation has been solved by HWCM with stability analysis [24]. Different types of nonlinear equations are solved using the Kudryashov method [25] and quintic B-splines collocation methods [26–29] including fractional Schrödinger equations [30], which have importance in applied mathematics and optics.The main objective of this study is to solve the Schrödinger equation via finite-difference formula. This Schrödinger equation is converted into system of algebraic equations using the difference formula at each steps. This technique is easy, reliable, and efficient to solve Equation (1) with optical soliton solutions. Due to the real and imaginary part of the complex function ux,t, Equation (1) is strenuous to tackle numerically, and most of the methods do not work due to the presence of imaginary part of the solution, such as Radial basis function methods.


## 2. Crank–Nicolson Scheme for Schrödinger Equations

Crank and Phyllis Nicolson (1947) proposed a method for the numerical solution of partial differential equations known as Crank–Nicolson method. The beauty of the method is the convergent and stability of results for all finite values ofλ, i.e., Δt/Δx2=λ [31].We also implement the Crank–Nicolson scheme to solve the time-dependent Schrödinger equation, which has a lot of applications in physics including optics and acoustics. The nondimensionilized form of Equation (1) can be written as(2)ι∂u∂t+η∂2u∂x2+αu2u+βfx,t=0,x∈a,b,t∈0,T,ι=−1,under initial condition,(3)ux,o=ϕx,and the boundary conditions,(4)ua,t=ξ0x,ub,t=ξ1x,where η≠0. If α=0, then Equation (2) becomes(5)ι∂u∂t+η∂2u∂x2+βfx,t=0,where the parameter β may be positive or negative. The complex-valued function u depends on x as well as on t and ϕ is a prescribed smooth function which represents the initial state. If η=−1 and β=0, then Equation (5) becomes(6)ι∂u∂t−∂2u∂x2=0.Using the finite-difference approximation, Equation (6) can be written as(7)ι∂u∂ti,j+1=12∂2u∂x2i,j+∂2u∂x2i,j+1ιUij+1−UijΔt=12Ui−1j−2Uij+Ui+1jΔx2+Ui−1j+1−2Uij+1+Ui+1j+1Δx2.Rearranging, we obtain(8)−λUi−1j+1+2ι+2λUij+1−λUi+1j+1=λUi−1j+2ι−2λUij+λUi+1j,where λ=Δt/Δx2, i=1,2,3,…,N and j=0,1,2,…,M. In Equation (7), we used the following initial and boundary conditions:(9)Ui0=ϕxi,i=1,2,3,…,N,U1j=ξ1tj,UNj=ξ2tj,j=0,1,2,…,M.In general, the left side of Equation (7) contains three unknown and right side contains three known values of U for each values of j.IntroducingN mesh points for each j such as i=1,2,3,…,N, then Equation (7) gives N−1 simultaneous equations for the N−1 unknown values for each time level.N−1For each j=1,2,3,…,M,(10)i=1:−λU01+2ι+2λU11−λU21=λU00+2ι−2λU10+λU20i=2:−λU11+2ι+2λU21−λU31=λU10+2ι−2λU20+λU30i=3:−λU21+2ι+2λU31−λU41=λU20+2ι−2λU30+λU40i=4:−λU312ι+2λ+U41−λU51=λU30+2ι−2λU40+λU50⋮⋮⋮i=N−1:−λUN−212ι+2λ+UN−11−λUN1=λUN−20+2ι−2λUN−10+λUN0.The above system of equation can be written in matrix form(11)AUj+1=BUj+Cj,where(12)A=2ι+2λ−λ00…0−λ2ι+2λ−λ0…00−λ2ι+2λ−λ…000−λ2ι+2λ−λ0⋮⋮⋮⋱⋱⋱0000−λ2ι+2λ,B=2ι−2λλ00…0λ2ι−2λλ0…00λ2ι−2λλ…∞000λ2ι−2λλ0⋮⋮⋮⋱⋱⋱0000λ2ι−2λ,Cj=λU0j+λU0j+10⋮⋮0λUNj+λUNj+1,Uj+1=U1j+1U2j+1⋮⋮UN−2j+1UN−1j+1,Uj=U1jU2j⋮⋮UN−2jUN−1j.Hence, we can solve Equation (8) for Uj+1:(13)Uj+1=A−1BUj+A−1Cj.Similarly, evaluation forj=2,3,…,M can also be made to get the required solution, where T=MΔt.


## 3. Result and Discussion

Various test problems have been explored in this section to demonstrate the accuracy and applicability of the method. The concerned numerical results are compared with the exact solution. For numerical simulation, we have used MATLAB R2013a software. To observe the performance of the current technique, we used the maximum absolute error normE∞, which is defined as(14)E∞=maxEabℐ,(15)Eabℐ=ℐexact−ℐN,ℐ=Reu+Imu,where exact and N represent the exact and numerical values, respectively, having real part Re and imaginary part Im.Test Problem 1.
Taking the following equation in the intervalx∈−2,2,(16)ι∂u∂t−∂2u∂x2=0.
The initial and boundary conditions are(17)ux,0=sinx,u−2,t=eιtsin−2,u2,t=eιtsin2.
The analytical solution isux,t=eιtsinx.
The numerical results are compared to the analytical solution in Figure1 for different values of collocation points, and it is obvious that increasing the number of collocation points improves the results. The findings are also shown in Figure 2 for various final times.
The numerical solutions are presented for different points inside the interval in Table1. In Table 2, the numerical results are displayed for various time iterations, and in the last row of the table, the exact solutions are also mentioned for comparison purposes. The absolute errors decrease as the number of meshes increases in Table 3. It can be observed from the table that raising the value of M has resulted in reasonably excellent accuracy.Figure 1
Comparison of exact and numerical solution atΔt=0.0013, Δx=0.1000, a=−2, b=2, and t=1 of test 1.Figure 2
Numerical result at differentt, Δt=0.0013, Δx=0.2, a=−2,b=2, and M=20 for test problem 1.Table 1
Numerical results for test problem1 at M=5, Δt=0.05, Δx=0.8, a=−2, and b=2.
ti=0i=1i=2i=3i=4i=5x=−2.0000x=−1.2000x=−0.4000x=0.4000x=1.2000x=2.0000j=00−0.9082−0.0454ι−0.9320−0.0000ι−0.3894−0.0000ι0.3894+0.0000ι0.9320+0.0000ι0.9082+0.0454ιj=10.0500−0.9048−0.0908ι−0.9363−0.0449ι−0.3890−0.0182ι0.3890+0.0182ι0.9363+0.0449ι0.9048+0.0908ιj=20.1000−0.8991−0.1359ι−0.9385−0.0908ι−0.3880−0.0361ι0.3880+0.0361ι0.9385+0.0908ι0.8991+0.1359ιj=30.1500−0.8912−0.1806ι−0.9383−0.1375ι−0.3864−0.0535ι0.3864+0.0535ι0.9383+0.1375ι0.8912+0.1806j=40.2000−0.8810−0.2250ι−0.9355−0.1849ι−0.3844−0.0706ι0.3844+0.0706ι0.9355+0.1849ι0.8810+0.2250ιj=50.2500−0.8687−0.2687ι−0.9300−0.2329ι−0.3822−0.0876ι0.3822+0.0876ι0.9300+0.2329ι0.8687+0.2687ιj=60.3000−0.8542−0.3118ι−0.9217−0.2810ι−0.3798−0.1046ι0.3798+0.1046ι0.9217+0.2810ι0.8542+0.3118ιj=70.3500−0.8375−0.3541ι−0.9106−0.3292ι−0.3771−0.1217ι0.3771+0.1217ι0.9106+0.3292ι0.8375+0.3541ι⋮⋮⋮⋮⋮⋮⋮j=100.4500−0.7752−0.4753ι−0.8597−0.4709ι−0.3667−0.1755ι0.3667+0.1755ι0.8597+0.4709ι0.7752+0.4753ι⋮⋮⋮⋮⋮⋮⋮j=140.7000−0.6653−0.6198ι−0.7544−0.6447ι−0.3402−0.2553ι0.3402+0.2553ι0.7544+0.6447ι0.6653+0.6198ι⋮⋮⋮⋮⋮⋮⋮j=180.9000−0.5289−0.7396ι−0.6141−0.7932ι−0.2875−0.3372ι0.2875+0.3372ι0.6141+0.7932ι0.5289+0.7396ι⋮⋮⋮⋮⋮⋮⋮j=201−0.4524−0.7887ι−0.5331−0.8568ι−0.2497−0.3736ι0.2497+0.3736ι0.5331+0.8568ι0.4524+0.7887ιExact solutiont=1−0.4913−0.7651ι−0.5036−0.7843ι−0.2104−0.3277ι0.2104+0.3277ι0.5036+0.7843ι0.4913+0.7651ιTable 2
Numerical values of the exact and approximate solution at domain−2,2, M=20, t=1, Δt=0.0013, and Δx=0.2 for test 1.
xuexactuN−2.0000−0.4913+0.7651ι−0.4903+0.7658ι−1.8000−0.5262+0.8195ι−0.5263+0.8218ι−1.6000−0.5401+0.8411ι−0.5413+0.8439ι−1.4000−0.5324+0.8292ι−0.5344+0.8324ι−1.2000−0.5036+0.7843ι−0.5058+0.7872ι−1.0000−0.4546+0.7081ι−0.4570+0.7105ι−0.8000−0.3876+0.6036ι−0.3900+0.6058ι−0.6000−0.3051+0.4751ι−0.3075+0.4769ι−0.4000−0.2104+0.3277ι−0.2126+0.3284ι−0.2000−0.1073+0.1672ι−0.1083+0.1671ι00.0000+0.0000ι0.0000+0.0000ι0.20000.1073−0.1672ι0.1083−0.1671ι0.40000.2104−0.3277ι0.2126−0.3284ι0.60000.3051−0.4751ι0.3075−0.4769ι0.80000.3876−0.6036ι0.3900−0.6058ι1.00000.4546−0.7081ι0.4570−0.7105ι1.20000.5036−0.7843ι0.5058−0.7872ι1.40000.5324−0.8292ι0.5344−0.8324ι1.60000.5401−0.8411ι0.5413−0.8439ι1.80000.5262−0.8195ι0.5263−0.8218ι2.00000.4913−0.7651ι0.4903−0.7658ιTable 3
Error at different values ofM at Δt=0.0013, a=−2, b=2, and t=1 for test 1.
MΔxE∞ReuE∞Imu50.80004.14×10−21.57×10−2100.40001.13×10−21.38×10−2150.26674.80×10−38.17×10−3200.20002.50×10−35.98×10−3250.16001.70×10−32.82×10−3300.13331.50×10−31.38×10−3350.11431.30×10−31.35×10−3400.10001.10×10−31.32×10−3Test problem 2.
Considering again the Schro..dinger equation,(18)ι∂u∂t−∂2u∂x2=0,x∈−π,π,where the required information can be retrieved from the analytical solution ux,t=eιtcosx. In Table 4, the maximum absolute errors are presented and the error decreases by increasing the collocation points M. Results are compared with the analytical solution in Figure 3 and also calculated for various time t in Figure 4. We also illustrated exact and numerical results with different input values in Figure 5.Table 4
Maximum error at different values ofM, Δt=0.0011, a=−π, b=π, and t=1 for test problem 2.
MΔxE∞ReuE∞Imu100.62833.09×10−24.68×10−2150.41891.35×10−22.52×10−2200.31427.90×10−31.02×10−2300.20943.80×10−38.28×10−3400.15711.80×10−33.84×10−3500.12571.20×10−32.62×10−3600.10471.20×10−31.60×10−3700.08981.13×10−31.34×10−3Figure 3
Numerical results atΔt=0.0011, Δx=0.3142, a=−π, b=π, and M=20 for test problem 1.Figure 4
Comparison of exact and numerical solution atΔt=0.0011, Δx=0.1571, a=−π, b=π, and t=1 of test problem 2.Figure 5
Comparison of the exact and numerical results atM=20, a=0, b=π, Δx=0.1571, and Δt=0.0014 for test problem 3.Test problem 3.
Finally, considering the following equation,(19)ι∂u∂t−∂2u∂x2=0,x∈0,π,with exact solution ux,t=eιt+x. In Figure 6, the numerical as well as the analytical solutions are presented in terms of both real and imaginary parts. It can be seen from this figure that the results easily approach to the exact solution for all the collocation points. Results are also calculated for various time t in Figure 6.Figure 6
Comparison of exact and numerical solution atΔt=0.0014, a=0, b=π, Δx=0.1571, and t=1 for test problem 3.


## 4. Conclusion

In this study, optical soliton solutions for the linear Schrödinger equation are conveyed by using the well-known Crank–Nicolson method. For this motivation, we use the Crank–Nicolson method for finding the numerical solution of Schrödinger equation. Three different examples are studied to verify the numerical results. From the comparison of the numerical solution with exact solution, it can be concluded that the numerical solutions converge to the exact solution without any restrictions ofΔt and Δx, but the accuracy increases by increasing the number of collocation point M; in other way, the results get accurate by decreasing Δt and Δx. The current approach can also be implemented to nonlinear Schrödinger equation and is the focus of our future task.



---

*Source: 6991067-2022-04-14.xml*