![]() |
Reference documentation for deal.II version 9.3.3
|
This tutorial depends on step-26, step-29.
This program was contributed by Wolfgang Bangerth (Colorado State University) and Yong-Yong Cai (Beijing Computational Science Research Center, CSRC) and is the result of the first author's time as a visitor at CSRC.
This material is based upon work partially supported by National Science Foundation grants OCI-1148116, OAC-1835673, DMS-1821210, and EAR-1925595; and by the Computational Infrastructure in Geodynamics initiative (CIG), through the National Science Foundation under Award No. EAR-1550901 and The University of California-Davis.
The Nonlinear Schrödinger Equation (NLSE) for a function and a potential
is a model often used in quantum mechanics and nonlinear optics. If one measures in appropriate quantities (so that
), then it reads as follows:
If there is no potential, i.e. , then it can be used to describe the propagation of light in optical fibers. If
, the equation is also sometimes called the Gross-Pitaevskii equation and can be used to model the time dependent behavior of Bose-Einstein condensates.
For this particular tutorial program, the physical interpretation of the equation is not of much concern to us. Rather, we want to use it as a model that allows us to explain two aspects:
At first glance, the equations appear to be parabolic and similar to the heat equation (see step-26) as there is only a single time derivative and two spatial derivatives. But this is misleading. Indeed, that this is not the correct interpretation is more easily seen if we assume for a moment that the potential and
. Then we have the equation
If we separate the solution into real and imaginary parts, , with
, then we can split the one equation into its real and imaginary parts in the same way as we did in step-29:
Not surprisingly, the factor in front of the time derivative couples the real and imaginary parts of the equation. If we want to understand this equation further, take the time derivative of one of the equations, say
(where we have assumed that, at least in some formal sense, we can commute the spatial and temporal derivatives), and then insert the other equation into it:
This equation is hyperbolic and similar in character to the wave equation. (This will also be obvious if you look at the video in the "Results" section of this program.) Furthermore, we could have arrived at the same equation for as well. Consequently, a better assumption for the NLSE is to think of it as a hyperbolic, wave-propagation equation than as a diffusion equation such as the heat equation. (You may wonder whether it is correct that the operator
appears with a positive sign whereas in the wave equation,
has a negative sign. This is indeed correct: After multiplying by a test function and integrating by parts, we want to come out with a positive (semi-)definite form. So, from
we obtain
. Likewise, after integrating by parts twice, we obtain from
the form
. In both cases do we get the desired positive sign.)
The real NLSE, of course, also has the terms and
. However, these are of lower order in the spatial derivatives, and while they are obviously important, they do not change the character of the equation.
In any case, the purpose of this discussion is to figure out what time stepping scheme might be appropriate for the equation. The conclusions is that, as a hyperbolic-kind of equation, we need to choose a time step that satisfies a CFL-type condition. If we were to use an explicit method (which we will not), we would have to investigate the eigenvalues of the matrix that corresponds to the spatial operator. If you followed the discussions of the video lectures (See also video lecture 26, video lecture 27, video lecture 28.) then you will remember that the pattern is that one needs to make sure that where
is the time step,
the mesh width, and
are the orders of temporal and spatial derivatives. Whether you take the original equation (
) or the reformulation for only the real or imaginary part, the outcome is that we would need to choose
if we were to use an explicit time stepping method. This is not feasible for the same reasons as in step-26 for the heat equation: It would yield impractically small time steps for even only modestly refined meshes. Rather, we have to use an implicit time stepping method and can then choose a more balanced
. Indeed, we will use the implicit Crank-Nicolson method as we have already done in step-23 before for the regular wave equation.
If one thought of the NLSE as an ordinary differential equation in which the right hand side happens to have spatial derivatives, i.e., write it as
one may be tempted to "formally solve" it by integrating both sides over a time interval and obtain
Of course, it's not that simple: the in the integrand is still changing over time in accordance with the differential equation, so we cannot just evaluate the integral (or approximate it easily via quadrature) because we don't know
. But we can write this with separate contributions as follows, and this will allow us to deal with different terms separately:
The way this equation can now be read is as follows: For each time interval , the change
in the solution consists of three contributions:
Operator splitting is now an approximation technique that allows us to treat each of these contributions separately. (If we want: In practice, we will treat the first two together, and the last one separate. But that is a detail, conceptually we could treat all of them differently.) To this end, let us introduce three separate "solutions":
These three "solutions" can be thought of as satisfying the following differential equations:
In other words, they are all trajectories that start at
and integrate up the effects of exactly one of the three terms. The increments resulting from each of these terms over our time interval are then
,
, and
.
It is now reasonable to assume (this is an approximation!) that the change due to all three of the effects in question is well approximated by the sum of the three separate increments:
This intuition is indeed correct, though the approximation is not exact: the difference between the exact left hand side and the term (i.e., the difference between the exact increment for the exact solution
when moving from
to
, and the increment composed of the three parts on the right hand side), is proportional to
. In other words, this approach introduces an error of size
. Nothing we have done so far has discretized anything in time or space, so the overall error is going to be
plus whatever error we commit when approximating the integrals (the temporal discretization error) plus whatever error we commit when approximating the spatial dependencies of
(the spatial error).
Before we continue with discussions about operator splitting, let us talk about why one would even want to go this way? The answer is simple: For some of the separate equations for the , we may have ways to solve them more efficiently than if we throw everything together and try to solve it at once. For example, and particularly pertinent in the current case: The equation for
, i.e.,
or equivalently,
can be solved exactly: the equation is solved by
This is easy to see if (i) you plug this solution into the differential equation, and (ii) realize that the magnitude is constant, i.e., the term
in the exponent is in fact equal to
. In other words, the solution of the ODE for
only changes its phase, but the magnitude of the complex-valued function
remains constant. This makes computing
particularly convenient: we don't actually need to solve any ODE, we can write the solution down by hand. Using the operator splitting approach, none of the methods to compute
therefore have to deal with the nonlinear term and all of the associated unpleasantries: we can get away with solving only linear problems, as long as we allow ourselves the luxury of using an operator splitting approach.
Secondly, one often uses operator splitting if the different physical effects described by the different terms have different time scales. Imagine, for example, a case where we really did have some sort of diffusion equation. Diffusion acts slowly, but if is large, then the "phase rotation" by the term
acts quickly. If we treated everything together, this would imply having to take rather small time steps. But with operator splitting, we can take large time steps
for the diffusion, and (assuming we didn't have an analytic solution) use an ODE solver with many small time steps to integrate the "phase rotation" equation for
from
to
. In other words, operator splitting allows us to decouple slow and fast time scales and treat them differently, with methods adjusted to each case.
While the method above allows to compute the three contributions in parallel, if we want, the method can be made slightly more accurate and easy to implement if we don't let the trajectories for the
start all at
, but instead let the trajectory for
start at the end point of the trajectory for
, namely
; similarly, we will start the trajectory for
start at the end point of the trajectory for
, namely
. This method is then called "Lie splitting" and has the same order of error as the method above, i.e., the splitting error is
.
This variation of operator splitting can be written as follows (carefully compare the initial conditions to the ones above):
(Obviously, while the formulas above imply that we should solve these problems in this particular order, it is equally valid to first solve for trajectory 3, then 2, then 1, or any other permutation.)
The integrated forms of these equations are then
From a practical perspective, this has the advantage that we need to keep around fewer solution vectors: Once has been computed, we don't need
any more; once
has been computed, we don't need
any more. And once
has been computed, we can just call it
because, if you insert the first into the second, and then into the third equation, you see that the right hand side of
now contains the contributions of all three physical effects:
(Compare this again with the "exact" computation of : It only differs in how we approximate
in each of the three integrals.) In other words, Lie splitting is a lot simpler to implement that the original method outlined above because data handling is so much simpler.
As mentioned above, Lie splitting is only accurate. This is acceptable if we were to use a first order time discretization, for example using the explicit or implicit Euler methods to solve the differential equations for
. This is because these time integration methods introduce an error proportional to
themselves, and so the splitting error is proportional to an error that we would introduce anyway, and does not diminish the overall convergence order.
But we typically want to use something higher order – say, a Crank-Nicolson or BDF2 method – since these are often not more expensive than a simple Euler method. It would be a shame if we were to use a time stepping method that is , but then lose the accuracy again through the operator splitting.
This is where the Strang splitting method comes in. It is easier to explain if we had only two parts, and so let us combine the effects of the Laplace operator and of the potential into one, and the phase rotation into a second effect. (Indeed, this is what we will do in the code since solving the equation with the Laplace equation with or without the potential costs the same – so we merge these two steps.) The Lie splitting method from above would then do the following: It computes solutions of the following two ODEs,
and then uses the approximation . In other words, we first make one full time step for physical effect one, then one full time step for physical effect two. The solution at the end of the time step is simply the sum of the increments due to each of these physical effects separately.
In contrast, Gil Strang (one of the titans of numerical analysis starting in the mid-20th century) figured out that it is more accurate to first do one half-step for one physical effect, then a full time step for the other physical effect, and then another half step for the first. Which one is which does not matter, but because it is so simple to do the phase rotation, we will use this effect for the half steps and then only need to do one spatial solve with the Laplace operator plus potential. This operator splitting method is now accurate. Written in formulas, this yields the following sequence of steps:
As before, the first and third step can be computed exactly for this particular equation, yielding
This is then how we are going to implement things in this program: In each time step, we execute three steps, namely
This structure will be reflected in an obvious way in the main time loop of the program.
From the discussion above, it should have become clear that the only partial differential equation we have to solve in each time step is
This equation is linear. Furthermore, we only have to solve it from to
, i.e., for exactly one time step.
To do this, we will apply the second order accurate Crank-Nicolson scheme that we have already used in some of the other time dependent codes (specifically: step-23 and step-26). It reads as follows:
Here, the "previous" solution (or the "initial
condition" for this part of the time step) is the output of the first phase rotation half-step; the output of the current step will be denoted by
.
is the length of the time step. (One could argue whether
and
live at time step
or
and what their upper indices should be. This is a philosophical discussion without practical impact, and one might think of
as something like
, and
as
if that helps clarify things – though, again
is not to be understood as "one third time step after
\_form#392" but more like "we've already done one third of the work necessary
for time step \_form#2675".)
If we multiply the whole equation with and sort terms with the unknown
to the left and those with the known
to the right, then we obtain the following (spatial) partial differential equation that needs to be solved in each time step:
As mentioned above, the previous tutorial program dealing with complex-valued solutions (namely, step-29) separated real and imaginary parts of the solution. It thus reduced everything to real arithmetic. In contrast, we here want to keep things complex-valued.
The first part of this is that we need to define the discretized solution as where the
are the usual shape functions (which are real valued) but the expansion coefficients
at time step
are now complex-valued. This is easily done in deal.II: We just have to use Vector<std::complex<double>> instead of Vector<double> to store these coefficients.
Of more interest is how to build and solve the linear system. Obviously, this will only be necessary for the second step of the Strang splitting discussed above, with the time discretization of the previous subsection. We obtain the fully discrete version through straightforward substitution of by
and multiplication by a test function:
or written in a more compact way:
Here, the matrices are defined in their obvious ways:
Note that all matrices individually are in fact symmetric, real-valued, and at least positive semidefinite, though the same is obviously not true for the system matrix and the corresponding matrix
on the right hand side.
The only remaining important question about the solution procedure is how to solve the complex-valued linear system
with the matrix and a right hand side that is easily computed as the product of a known matrix and the previous part-step's solution. As usual, this comes down to the question of what properties the matrix
has. If it is symmetric and positive definite, then we can for example use the Conjugate Gradient method.
Unfortunately, the matrix's only useful property is that it is complex symmetric, i.e., , as is easy to see by recalling that
are all symmetric. It is not, however, Hermitian, which would require that
where the bar indicates complex conjugation.
Complex symmetry can be exploited for iterative solvers as a quick literature search indicates. We will here not try to become too sophisticated (and indeed leave this to the Possibilities for extensions section below) and instead simply go with the good old standby for problems without properties: A direct solver. That's not optimal, especially for large problems, but it shall suffice for the purposes of a tutorial program. Fortunately, the SparseDirectUMFPACK class allows solving complex-valued problems.
Initial conditions for the NLSE are typically chosen to represent particular physical situations. This is beyond the scope of this program, but suffice it to say that these initial conditions are (i) often superpositions of the wave functions of particles located at different points, and that (ii) because corresponds to a particle density function, the integral
corresponds to the number of particles in the system. (Clearly, if one were to be physically correct, better be a constant if the system is closed, or
if one has absorbing boundary conditions.) The important point is that one should choose initial conditions so that
makes sense.
What we will use here, primarily because it makes for good graphics, is the following:
where is the distance from the (fixed) locations
, and
are chosen so that each of the Gaussians that we are adding up adds an integer number of particles to
. We achieve this by making sure that
is a positive integer. In other words, we need to choose as an integer multiple of
assuming for the moment that – which is of course not the case, but we'll ignore the small difference in integral.
Thus, we choose for all, and
. This
is small enough that the difference between the exact (infinite) integral and the integral over
should not be too concerning. We choose the four points
as
– also far enough away from the boundary of
to keep ourselves on the safe side.
For simplicity, we pose the problem on the square . For boundary conditions, we will use time-independent Neumann conditions of the form
This is not a realistic choice of boundary conditions but sufficient for what we want to demonstrate here. We will comment further on this in the Possibilities for extensions section below.
Finally, we choose , and the potential as
Using a large potential makes sure that the wave function remains small outside the circle of radius 0.7. All of the Gaussians that make up the initial conditions are within this circle, and the solution will mostly oscillate within it, with a small amount of energy radiating into the outside. The use of a large potential also makes sure that the nonphysical boundary condition does not have too large an effect.
The program starts with the usual include files, all of which you should have seen before by now:
Then the usual placing of all content of this program into a namespace and the importation of the deal.II namespace into the one we will work in:
NonlinearSchroedingerEquation
classThen the main class. It looks very much like the corresponding classes in step-4 or step-6, with the only exception that the matrices and vectors and everything else related to the linear system are now storing elements of type std::complex<double>
instead of just double
.
Before we go on filling in the details of the main class, let us define the equation data corresponding to the problem, i.e. initial values, as well as a right hand side class. (We will reuse the initial conditions also for the boundary values, which we simply keep constant.) We do so using classes derived from the Function class template that has been used many times before, so the following should not look surprising. The only point of interest is that we here have a complex-valued problem, so we have to provide the second template argument of the Function class (which would otherwise default to double
). Furthermore, the return type of the value()
functions is then of course also complex.
What precisely these functions return has been discussed at the end of the Introduction section.
NonlinearSchroedingerEquation
classWe start by specifying the implementation of the constructor of the class. There is nothing of surprise to see here except perhaps that we choose quadratic ( ) Lagrange elements – the solution is expected to be smooth, so we choose a higher polynomial degree than the bare minimum.
The next function is the one that sets up the mesh, DoFHandler, and matrices and vectors at the beginning of the program, i.e. before the first time step. The first few lines are pretty much standard if you've read through the tutorial programs at least up to step-6:
Next, we assemble the relevant matrices. The way we have written the Crank-Nicolson discretization of the spatial step of the Strang splitting (i.e., the second of the three partial steps in each time step), we were led to the linear system . In other words, there are two matrices in play here – one for the left and one for the right hand side. We build these matrices separately. (One could avoid building the right hand side matrix and instead just form the action of the matrix on
in each time step. This may or may not be more efficient, but efficiency is not foremost on our minds for this program.)
Having set up all data structures above, we are now in a position to implement the partial steps that form the Strang splitting scheme. We start with the half-step to advance the phase, and that is used as the first and last part of each time step.
To this end, recall that for the first half step, we needed to compute . Here,
and
are functions of space and correspond to the output of the previous complete time step and the result of the first of the three part steps, respectively. A corresponding solution must be computed for the third of the part steps, i.e.
, where
is the result of the time step as a whole, and its input
is the result of the spatial step of the Strang splitting.
An important realization is that while may be a finite element function (i.e., is piecewise polynomial), this may not necessarily be the case for the "rotated" function in which we have updated the phase using the exponential factor (recall that the amplitude of that function remains constant as part of that step). In other words, we could compute
at every point
, but we can't represent it on a mesh because it is not a piecewise polynomial function. The best we can do in a discrete setting is to compute a projection or interpolation. In other words, we can compute
where
is a projection or interpolation operator. The situation is particularly simple if we choose the interpolation: Then, all we need to compute is the value of the right hand side at the node points and use these as nodal values for the vector
of degrees of freedom. This is easily done because evaluating the right hand side at node points for a Lagrange finite element as used here requires us to only look at a single (complex-valued) entry of the node vector. In other words, what we need to do is to compute
where
loops over all of the entries of our solution vector. This is what the function below does – in fact, it doesn't even use separate vectors for
and
, but just updates the same vector as appropriate.
The next step is to solve for the linear system in each time step, i.e., the second half step of the Strang splitting we use. Recall that it had the form where
and
are the matrices we assembled earlier.
The way we solve this here is using a direct solver. We first form the right hand side using the SparseMatrix::vmult() function and put the result into the
system_rhs
variable. We then call SparseDirectUMFPACK::solver() which takes as argument the matrix and the right hand side vector and returns the solution in the same vector
system_rhs
. The final step is then to put the solution so computed back into the solution
variable.
The last of the helper functions and classes we ought to discuss are the ones that create graphical output. The result of running the half and full steps for the local and spatial parts of the Strang splitting is that we have updated the solution
vector to the correct value at the end of each time step. Its entries contain complex numbers for the solution at the nodes of the finite element mesh.
Complex numbers are not easily visualized. We can output their real and imaginary parts, i.e., the fields and
, and that is exactly what the DataOut class does when one attaches as complex-valued vector via DataOut::add_data_vector() and then calls DataOut::build_patches(). That is indeed what we do below.
But oftentimes we are not particularly interested in real and imaginary parts of the solution vector, but instead in derived quantities such as the magnitude and phase angle
of the solution. In the context of quantum systems such as here, the magnitude itself is not so interesting, but instead it is the "amplitude",
that is a physical property: it corresponds to the probability density of finding a particle in a particular place of state. The way to put computed quantities into output files for visualization – as used in numerous previous tutorial programs – is to use the facilities of the DataPostprocessor and derived classes. Specifically, both the amplitude of a complex number and its phase angles are scalar quantities, and so the DataPostprocessorScalar class is the right tool to base what we want to do on.
Consequently, what we do here is to implement two classes ComplexAmplitude
and ComplexPhase
that compute for each point at which DataOut decides to generate output, the amplitudes and phases
of the solution for visualization. There is a fair amount of boiler-plate code below, with the only interesting parts of the first of these two classes being how its
evaluate_vector_field()
function computes the computed_quantities
object.
(There is also the rather awkward fact that the std::norm() function does not compute what one would naively imagine, namely , but returns
instead. It's certainly quite confusing to have a standard function mis-named in such a way...)
The second of these postprocessor classes computes the phase angle of the complex-valued solution at each point. In other words, if we represent , then this class computes
. The function std::arg does this for us, and returns the angle as a real number between
and
.
For reasons that we will explain in detail in the results section, we do not actually output this value at each location where output is generated. Rather, we take the maximum over all evaluation points of the phase and then fill each evaluation point's output field with this maximum – in essence, we output the phase angle as a piecewise constant field, where each cell has its own constant value. The reasons for this will become clear once you read through the discussion further down below.
Having so implemented these post-processors, we create output as we always do. As in many other time-dependent tutorial programs, we attach flags to DataOut that indicate the number of the time step and the current simulation time.
The remaining step is how we set up the overall logic for this program. It's really relatively simple: Set up the data structures; interpolate the initial conditions onto finite element space; then iterate over all time steps, and on each time step perform the three parts of the Strang splitting method. Every tenth time step, we generate graphical output. That's it.
The rest is again boiler plate and exactly as in almost all of the previous tutorial programs:
Running the code results in screen output like the following:
Running the program also yields a good number of output files that we will visualize in the following.
The output_results()
function of this program generates output files that consist of a number of variables: The solution (split into its real and imaginary parts), the amplitude, and the phase. If we visualize these four fields, we get images like the following after a few time steps (at time , to be precise:
While the real and imaginary parts of the solution shown above are not particularly interesting (because, from a physical perspective, the global offset of the phase and therefore the balance between real and imaginary components, is meaningless), it is much more interesting to visualize the amplitude and phase
of the solution and, in particular, their evolution. This leads to pictures like the following:
The phase picture shown here clearly has some flaws:
nic_Edge
map in which both of the extreme values are shown as black.ComplexPhase
class just uses the maximal phase angle encountered on each cell.With these modifications, the phase plot now looks as follows:
Finally, we can generate a movie out of this. (To be precise, the video uses two more global refinement cycles and a time step half the size of what is used in the program above.) The author of these lines made the movie with VisIt, because that's what he's more familiar with, and using a hacked color map that is also cyclic – though this color map lacks all of the skill employed by the people who wrote the posts mentioned in the links above. It does, however, show the character of the solution as a wave equation if you look at the shaded part of the domain outside the circle of radius 0.7 in which the potential is zero – you can see how every time one of the bumps (showing the amplitude ) bumps into the area where the potential is large: a wave travels outbound from there. Take a look at the video:
So why did I end up shading the area where the potential is large? In that outside region, the solution is relatively small. It is also relatively smooth. As a consequence, to some approximate degree, the equation in that region simplifies to
or maybe easier to read:
To the degree to which this approximation is valid (which, among other things, eliminates the traveling waves you can see in the video), this equation has a solution
Because is large, this means that the phase rotates quite rapidly. If you focus on the semi-transparent outer part of the domain, you can see that. If one colors this region in the same way as the inner part of the domain, this rapidly flashing outer part may be psychedelic, but is also distracting of what's happening on the inside; it's also quite hard to actually see the radiating waves that are easy to see at the beginning of the video.
The solver chosen here is just too simple. It is also not efficient. What we do here is give the matrix to a sparse direct solver in every time step and let it find the solution of the linear system. But we know that we could do far better:
In order to be usable for actual, realistic problems, solvers for the nonlinear Schrödinger equation need to utilize boundary conditions that make sense for the problem at hand. We have here restricted ourselves to simple Neumann boundary conditions – but these do not actually make sense for the problem. Indeed, the equations are generally posed on an infinite domain. But, since we can't compute on infinite domains, we need to truncate it somewhere and instead pose boundary conditions that make sense for this artificially small domain. The approach widely used is to use the Perfectly Matched Layer method that corresponds to a particular kind of attenuation. It is, in a different context, also used in step-62.
Finally, we know from experience and many other tutorial programs that it is worthwhile to use adaptively refined meshes, rather than the uniform meshes used here. It would, in fact, not be very difficult to add this here: It just requires periodic remeshing and transfer of the solution from one mesh to the next. step-26 will be a good guide for how this could be implemented.