A Mesh-free Numerical Method for three-dimensional Nonlinear Schrödinger Equation

1. A Mesh-free Numerical Method for three-dimensionalNonlinear Schrödinger Equation Department of Computer Science and Information Systems Birkbeck, University of London Thomas C.L. Yue tclyue@gmail.com Feb 09, 2011

2. Overview • Physical motivation of the problem • Dimensionless Gross-Pitaevskii equation (GPE) • Introduction to Radial basis functions (RBF) • Global supported strictly positive definite radial basis functions • Compactly supported radial basis funtions • Kansa’s method (asymmetric collocation) • Meshfree solution of cubic Nonlinear Schrodinger Equation • Numerical experiments and validation

3. Physical Motivation

4. Physical Motivation History of Bose Einstein Condensation (BEC) [1,2] • First predicted by Bose & Einstein (1924) • Experimentally observed in University of Colorado JILA lab (1995) What is BEC? [1,2] • A phase of matter where all particles occupy the same quantum state • Occurs when diulated bosons (integer spin particles) gas are cooled to extremely low temperature (10-9K) • Individual particle wave functions behave as a single wave function

5. Physical Motivation 1. High temperature particle behaviour dominated 2. Low temperature λdB α T -0.5 3. T=Tcrit Bose Einstein Condensate 4. T=0 Giant Matter Wave Fig1.A visual description of how a gas of bosonic-atoms behave at various temperatures (T). [1]

6. Experimental Results of BEC JILA (95’,Rb,5,000) ETH (02’,Rb, 300,000)

7. Gross–Pitaevskii equation • Hartree–Fock approximation [1,2] • The many-body wavefunction is written as productsof individual wave functions of each bosons [1,2] • The Hamiltonian • The conserved quantities

8. Gross–Pitaevskii equation • At temperature T<<Tcirt the dynamics of BEC is modeled the Gross–Pitaevskii equation [1,2] • Dimensionless variables introduced by Bao et al. (2003) [3]

9. Gross–Pitaevskii equation • Rearranging the equation and defining the following constants • The dimensionless Gross–Pitaevskii equation Note: This is mathematical equivalent to the cubic Nonlinear Schrödinger Equation (NLS)

10. Existing numerical methods for Nonlinear Schrödinger Equation Existing numerical methods for NLS Spectral Methods • Pseudo-spectral method (Muruganandam et al) • Time splitting Fourier spectral approximation (Bao et al.) • Split-step Fourier spectral method (Weideman) Mesh-based Methods • Galerkin spectral (Dion et al.) • Finite Element (Carl Joachim, Berdal Haga) • Split-step finite difference method (Wang)

11. Existing numerical methods for Nonlinear Schrödinger Equation Existing numerical methods for NLS Spectral Methods • Pseudo-spectral method (Muruganandam et al) • Time splitting Fourier spectral approximation (Bao et al.) • Split-step Fourier spectral method (Weideman) Mesh-based Methods • Galerkin spectral (Dion et al.) • Finite Element (Carl Joachim, Berdal Haga) • Split-step finite difference method (Wang) Require mesh generation and re-meshing

12. Radial Basis Functions

13. Radial basis function • What is a radial basis function (RBF)? [4,5]

14. RBF scattered data approximation • Given a set of data {x1...xN} and the corresponding known values {f(x1)..f(xN)}. Find the function f(x) that describes the data set. • Is the system guaranteed to be solvable? • Are the solutions unique?

15. RBF scattered data approximation Fig 2. Interpolation of f(x,y) with Gaussian RBF with c=1/3 and N=25. (left) shows the random generated data points, (mid) shows the centred at the collocation points, (right) shows the interpolated surface.

16. Background of radial basis functions • The system is solvable and unique provided the coefficient matrix is positive definite. [4,5,11]

17. Background of radial basis functions Globally supported strictly positive definite radial basis functions (GSRBF) • Leads to dense coefficient matrix • In many cases the coefficient matrix is ill-conditioned • For matrix inversion Schaback (2007) suggested • Singular Value Decomposition • Regularization techniques

18. Background of radial basis functions Compactly supported radial basis functions (CSRBF) • Wu and Wendland introduced the compactly supported RBF (CSRBF) [4,5] • Leads to sparse coefficient matrix • Reduce ill-conditioning of the resultant coefficient matrix • The usage of CSRBF will be explored in 3D NLS numerical experiment

19. Error Behaviour of RBF techniques • Trade off principle Schaback (1995) [5] Theorem: It is impossible to construct radial basis functions which guarantees good stability and small errors at the same time. • Driscoll and Fornberg (2002) observed the "Flat Limit” [6] c->∞ leads to highly ill-conditioned RBF interpolation matrix c->0 implies highly localized RBFs such that it fails to approximate data between collocation points

20. Error Behaviour of RBF techniques • Wright, Fornberg, Larsson (2004) [7] • With increasing shape parameter, interpolation error decreases sharply until the minimum numerical error is reached. • For any increasing shape parameter, interpolation error rapidly increases. The rapid decrease of interpolation error reaches a minimum.

21. Solving PDE with radial basis functions • Kansa (1990) proposed a direct approach to approximate the solution of PDE by • where Ф represents any RBF and p(x) is basis polynomial of up to order m. • Consider a linear PDE boundary value problem • where the linear operator L operates on the interior points Ω/∂Ω, the operator B specifies the boundary conditions for collocations on the boundaries ∂Ω.

22. Solving PDE with radial basis functions • Applying the RBF approximation the domain with Ni interior points in Ω/∂Ω and Nb boundary points on ∂Ω yields N equations • To remove the extra m degrees of freedom of the polynomial p(x)

23. Solving PDE with radial basis functions • Rewriting in matrix form • Note: The resultant PDE matrix is asymmetric. Hence Kansa method is also known as asymmetric collocation method.

24. Solving time-dependent PDE with θ-method and RBF • Some common methods for time-dependent PDE • θ-method • Runge-Kutta • Laplace Transform • θ-method • Based on the discretization of time-domain of the PDE. • The forward and backward time-step is weighted by (0≤θ≤1) • Consider the following time-dependent linear PDE problem

25. Solving time-dependent PDE with θ-method and RBF • constructing a time-domain mesh for M units, such that each time increment is denoted by tn=ndt, n=1..M, dt=T/M. • Hence the approximated PDE problem becomes • Approximate spatial variables by radial basis functions (ie. Kansa method)

26. Meshfree Numerical Method for Nonlinear Schrödinger Equation

27. Mesh-free Numerical Method for Nonlinear Schrödinger Equation • Recall: The equation for modelling dynamics of Bose-Einstein condensate (time-dependent Gross–Pitaevskii equation) • The Gross–Pitaevskii equation is mathematical equivalent to the cubic Nonlinear Schrödinger equation. • The parameter q controls the interaction between particles • q>0 defocusing interaction • q<0 focusing interaction

28. Mesh-free Numerical Method for Nonlinear Schrödinger Equation • The full 3D cubic Nonlinear Schrodinger equation (NLS) with initial and boundary conditions

29. Mesh-free Numerical Method for Nonlinear Schrödinger Equation • Key-steps for deriving the mesh-free method for NLS • separate the original NLS into real r(x,t) and imaginary parts s(x,t) • apply θ-method in time-domain • linearize PDE using the approach in Dereli (2009) • apply Kansa asymmetric collocation to spatial variables • Advantages of the proposed mathematical method • entirely meshfree • solves NLS in various dimensions d≤3 • flexible for selecting radial basis functions • easy to implement (~200 lines of matlab code)

30. Derivation of the proposed method • Separating the original NLS with respect to real r(x,t) and imaginary parts s(x,t) yields a system of PDEs. • Applying θ-method in time-domain

31. Derivation of the proposed method • Using the approach by Dereliet al (2009) [8] the variables (r*,s*) are introduced to approximate the solutions sufficient close to (rn+1,sn+1)

32. Derivation of the proposed method • Defining an auxiliary variable • Rewrite the real and imaginary parts of NLS using the definition of (r*,s*) and α: (Real) (Imaginary)

33. Derivation of the proposed method • Apply the RBF approximation to the real part r(x,t) and imaginary part s(x,t) of the wavefunction Ψ (x,t) and its spatial derivatives

34. Derivation of the proposed method

35. Derivation of the proposed method

36. Derivation of the proposed method • Final matrix form results a system of 2Nx2N equations • Solved via Singular Value Decomposition at each time-step to find RBF coefficients ζn+1 • Specific cases of θ-method • θ=0 explicit method • θ=0.5 semi-implicit method • θ=1 implicit method

37. Implementation flow-chart Set up physical geometries and potential function Compute initial conditions start Kernel of the method Assemble matrices for computation Visualize results while t<T start Update coefficients Conduct matrix inversion (compute new coefficients) Output numerical solution if(t==T)

38. Numerical Experiments

39. Radial basis functions in this project Globally supported strictly positive definite radial basis function (GSRBF) Compactly supported radial basis function (CSRBF)for 3D problem

40. 1D NLS numerical example • We consider a 1D test case in Deconinck et al. (2001) to model the stability of Bose Einstein Condensates and Wang (2005). [11]

41. 1D NLS numerical example • Comparison of absolute error between split-step finite difference method (SSFD) in Weideman (1986) and split-step Fourier spectral (SSFS) in Wang (2005). [11] Table 1. Absolute error comparison of RBF-θ and earlier methods. The solution is computed using RBF= Gaussian, θ=0.5, M=200, N=128, c=2.5. Table 2. Maximum relative error and maximum RMS error of real and imaginary parts of the wavefunction at T=1 generated by different globally supported strictly positive definite RBFs with M=500, N=128.

42. Fig 6. Real and imaginary parts of the numerical solution and the corresponding relative error at T=1 computed by RBF=Gaussian, M=500, N=128, c=2.5, θ=0.5.

43. Fig 7. Particle density (top) and relative error (bottom) of numerical solution at T=1 with M=500, N=128, c=2.5, θ=0.5, RBF=Gaussian.

44. 2D NLS numerical experiment • Consider a 2D defocusing interaction where q=1, k=1

45. 2D NLS numerical results Table 5. Maximum relative error, RMS error for different GSRBFs with M=2000, N=100, T=1. Table 6. Maximum relative error and RMS error of particle density at T=1 generated by different GSRBFs with M=2000, N=100.

46. Fig 10. Real and imaginary parts of numerical solutions and the corresponding relative error at time T=1 computed by M=2000, N=100, c=0.7, θ=1, RBF=Gaussian

47. Fig 11. Particle density (top) and relative error (bottom) of numerical solution at T=1 computed by M=2000, N=100, c=0.7, θ=1, RBF=Gaussian

48. 3D NLS numerical experiment • Consider a 3D focusing example where q=-1, k=2

49. 3D NLS numerical results Numerical results for all θ-methods and GSRBF combinations Table 7. Maximum relative error and RMS error of particle density at T=1 generated by various GSRBFs.

50. Fig 12. Real and imaginary parts of numerical solutions and the corresponding relative error at time T=1 computed by M=800, N=216, c=2.0, θ=1,RBF=IMQ.