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.
Introduction
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:
It is a complex-valued equation for . We have previously seen complex-valued equations in step-29, but there have opted to split the equations into real and imaginary parts and consequently ended up solving a system of two real-valued equations. In contrast, the goal here is to show how to solve problems in which we keep everything as complex numbers.
The equation is a nice model problem to explain how operator splitting methods work. This is because it has terms with fundamentally different character: on the one hand, is a regular spatial operator in the way we have seen many times before; on the other hand, has no spatial or temporal derivatives, i.e., it is a purely local operator. It turns out that we have efficient methods for each of these terms (in particular, we have analytic solutions for the latter), and that we may be better off treating these terms differently and separately. We will explain this in more detail below.
A note about the character of the equations
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.
The general idea of operator splitting
Note
The material presented here is also discussed in video lecture 30.25. (All video lectures are also available here.)
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:
The contribution of the Laplace operator.
The contribution of the potential .
The contribution of the "phase" term .
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.
Operator splitting: the "Lie splitting" approach
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.
Operator splitting: the "Strang splitting" approach
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
Update the solution value at each node by analytically integrating the phase rotation equation by one half time step;
Solving the space-time equation that corresponds to the full step for , namely , with initial conditions equal to the solution of the first half step above.
Update the solution value at each node by analytically integrating the phase rotation equation by another half time step.
This structure will be reflected in an obvious way in the main time loop of the program.
Time discretization
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#440" but more like "we've already done one third of the work necessary
for time step \_form#3124".)
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:
Spatial discretization and dealing with complex variables
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.
Linear solvers
Note
The material presented here is also discussed in video lecture 34. (All video lectures are also available here.)
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.
Definition of the test case
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 commented program
Include files
The program starts with the usual include files, all of which you should have seen before by now:
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/error_estimator.h>
#include <fstream>
#include <iostream>
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:
Then 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.
template <int dim>
class InitialValues : publicFunction<dim, std::complex<double>>
Implementation of the NonlinearSchroedingerEquation class
We 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.
Setting up data structures and assembling matrices
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.
Running the code results in screen output like the following:
Number of active cells: 4096
Number of degrees of freedom: 16641
Time step 1 at t=0
Time step 2 at t=0.00390625
Time step 3 at t=0.0078125
Time step 4 at t=0.0117188
[...]
Running the program also yields a good number of output files that we will visualize in the following.
Visualizing the solution
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:
First, phase is a "cyclic quantity", but the color scale uses a fundamentally different color for values close to than for values close to . This is a nuisance – what we need is a "cyclic color map" that uses the same colors for the two extremes of the range of the phase. Such color maps exist, see this blog post of Nicolás Guarín-Zapata or this StackExchange post, for example. The problem is that the author's favorite one of the two big visualization packages, VisIt, does not have any of these color maps built in. In an act of desperation, I therefore had to resort to using Paraview given that it has several of the color maps mentioned in the post above implemented. The picture below uses the nic_Edge map in which both of the extreme values are shown as black.
There is a problem on cells in which the phase wraps around. If at some evaluation point of the cell the phase value is close to and at another evaluation point it is close to , then what we would really like to happen is for the entire cell to have a color close to the extremes. But, instead, visualization programs produce a linear interpolation in which the values within the cell, i.e., between the evaluation points, is linearly interpolated between these two values, covering essentially the entire range of possible phase values and, consequently, cycling through the entire rainbow of colors from dark red to dark green over the course of one cell. The solution to this problem is to just output the phase value on each cell as a piecewise constant. Because averaging values close to the and is going to result in an average that has nothing to do with the actual phase angle, the 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.
Possibilities for extensions
Better linear solvers
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:
First, we should make use of the fact that the matrix doesn't actually change from time step to time step. This is an artifact of the fact that we here have constant boundary values and that we don't change the time step size – two assumptions that might not be true in actual applications. But at least in cases where this does happen to be the case, it would make sense to only factorize the matrix once (i.e., compute and factors once) and then use these factors for all following time steps until the matrix changes and requires a new factorization. The interface of the SparseDirectUMFPACK class allows for this.
Ultimately, however, sparse direct solvers are only efficient for relatively small problems, say up to a few 100,000 unknowns. Beyond this, one needs iterative solvers such as the Conjugate Gradient method (for symmetric and positive definite problems) or GMRES. We have used many of these in other tutorial programs. In all cases, they need to be accompanied by good preconditioners. For the current case, one could in principle use GMRES – a method that does not require any specific properties of the matrix – as the outer solver but at least at the time of writing this sentence (in 2022), the SolverGMRES class can only handle real-valued linear systems. This can be overcome by implementing a variation of GMRES that can deal with complex-valued matrices and vectors, see for example [Fraysse2005] . Even better would be to implement an iterative scheme that exploits the one structural feature we know is true for this problem: That the matrix is complex-symmetric (albeit not Hermitian), for which a literature search would probably find schemes as well.
A different strategy towards iterative solvers would be to break the linear system into a block system of real and imaginary components, like we did in step-29. This would then enable using real-valued iterative solvers on the outer level (e.g., the existing GMRES implementation), but one would have to come up with preconditioners that exploit the block structure. There is, again, literature on the topic, of which we simply point out a non-representative sample: [Axelsson2014] , [Day2001] , [Liao2016] .
Boundary conditions
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.
Adaptive meshes
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.