# A Reduced Order technique to study bifurcating phenomena: application to the Gross-Pitaevskii equation

We are interested in the steady state solution of the Gross-Pitaevskii equation, which describes the ground state of a quantum system of identical bosons under simplifying hypotheses, and their stability properties. Often referred to as a non-linear Schrödinger equation, the Gross-Pitaevskii equation models certain classes of Bose-Einstein condensates (BECs), a special state of matter formed by bosons at ultra-low temperatures. It is well known that the solutions of the Gross-Pitaevskii equation with a parabolic trap in two dimensions exhibit a rich bifurcating behavior, which includes symmetry-breaking bifurcations and vortex-bearing states. The bifurcation diagram plots the number of bosons in a BEC as a function of the chemical potential. In order to trace a bifurcation diagram, one has to compute multiple solution of a parametrized, non-linear problem. We will start with a one parameter study (the chemical potential) and then consider a two parameter case (chemical potential and the normalized trap strength). In this work, we propose to apply a Reduced Order Modeling (ROM) technique to reduce the demanding computational costs associated with this stability analysis.

## Authors

• 5 publications
• 18 publications
• 66 publications
12/06/2020

### Dynamics of solutions in the generalized Benjamin-Ono equation: a numerical study

We consider the generalized Benjamin-Ono (gBO) equation on the real line...
02/21/2021

### Bifurcation analysis of two-dimensional Rayleigh–Bénard convection using deflation

We perform a bifurcation analysis of the steady state solutions of Rayle...
02/27/2018

### Finding steady-state solutions for ODE systems of zero, first and homogeneous second-order chemical reactions is NP-hard

In the context of modeling of cell signaling pathways, a relevant step i...
12/30/2020

### Numerical study of soliton stability, resolution and interactions in the 3D Zakharov-Kuznetsov equation

We present a detailed numerical study of solutions to the Zakharov-Kuzne...
05/06/2021

### On Linear Damping Around Inhomogeneous Stationary States of the Vlasov-HMF Model

We study the dynamics of perturbations around an inhomogeneous stationar...
06/29/2021

### Tensor-train approximation of the chemical master equation and its application for parameter inference

In this work, we perform Bayesian inference tasks for the chemical maste...
05/01/2021

### Asymptotic behavior of fronts and pulses of the bidomain model

The bidomain model is the standard model for cardiac electrophysiology. ...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

We are interested in the steady state solution of the Gross–Pitaevskii equation, which describes the ground state of a quantum system of identical bosons under simplifying hypotheses, and their stability properties. Often referred to as a non-linear Schrödinger equation, the Gross–Pitaevskii equation models certain classes of Bose-Einstein condensates (BECs), a special state of matter formed by bosons at ultra-low temperatures. It is well known that the solutions of the Gross–Pitaevskii equation with a parabolic trap in two dimensions exhibit a rich bifurcating behavior [16, 18, 6], which includes symmetry-breaking bifurcations and vortex-bearing states [10].

The bifurcation diagram plots the number of bosons in a BEC as a function of the chemical potential. In order to trace a bifurcation diagram, one has to compute multiple solution of a parametrized, non-linear problem. We will start with a one parameter study (the chemical potential) and then consider a two parameter case (chemical potential and the normalized trap strength). Moreover, we are interested in studying the stability property of each steady solution to the Gross–Pitaevskii equation. We assume a perturbation around the steady state characterized by a “small” perturbation amplitude and a complex eigenfrequency. At order

, we obtain an eigenvalue problem that allows us to characterize the stability of the steady solution and, in the case of an unstable state, to classify the instability. Tracing the bifurcation diagram and the numerical study of the stability of the solutions can be every expensive in terms of computational time. In this work, we propose to apply a Reduced Order Modeling (ROM) technique to reduce the demanding computational costs associated with this stability analysis.

## 2 The Gross–Pitaevskii equation

The Gross–Pitaevskii equation models certain classes of Bose-Einstein condensates. A Bose-Einstein condensate (BEC) is a special state of matter formed by an unlimited number of bosons that “condense” into the same energy state at low temperatures. A BEC is formed by cooling a gas of extremely low density, about one-hundred-thousandth the density of normal air, to ultra-low temperatures (close to absolute zero).

A quantum system is the environment to be studied in terms of wave-particle duality (i.e., all particles exhibit a wave nature and viceversa). A quantum system involves the wave-function and its constituents, such as the momentum and wavelength. The Gross–Pitaevskii equation describes the ground state of a quantum system of identical bosons using two simplifications: the Hartree–Fock approximation and the pseudopotential interaction model. In the Hartree–Fock approximation, the total wave-function of a system of bosons is taken as a product of single-particle functions :

 Φtot(r1,r2,…,rN)=N∏i=1Φ(ri),

where is the coordinate of the -th boson. If the single-particle wave-function satisfies the Gross–Pitaevskii equation, the total wave-function minimizes the expectation value (i.e., the probabilistic expected value of the result of an experiment) of the pseudopotential model Hamiltonian under normalization condition:

 N=∫Dρ dr,ρ=|Φ|2, (1)

where is the domain under consideration and is interpreted as the particle density. The Gross–Pitaevskii equation reads: Find the single-particle wave-function such that

 i∂tΦ= −12ΔΦ+|Φ|2Φ+W(r)Φ% in D, (2)

where is the imaginary unit, is the radial coordinate, and is the external potential, with being the normalized trap strength, i.e. the ratio of trappings along and transverse to the plane. Here, we choose unless specified otherwise. Notice that we consider a single well potential. Eq. (2) is similar in form to the Ginzburg–Landau equation and is sometimes referred to as a non-linear Schrödinger equation. Obviously, eq. (2) needs to be supplemented with suitable boundary conditions.

The construction of the steady solution is based on the ansatz:

 Φ(r,t)=ϕ(r)exp(−iμt),ϕ(r):¯¯¯¯¯D→C, (3)

where is the chemical potential, which has to satisfy . By plugging (3) into (2), we obtain non-linear problem

 G(ϕ;μ)≐ −12Δϕ+|ϕ|2ϕ+W(r)ϕ−μϕ=0. (4)

It is well known that the solutions of one-dimensional version of problem (4) exhibit a bifurcating behavior [14, 13, 2, 8], which is not particularly rich though. The bifurcations occurring in the two-dimensional problem (4) are far more interesting [16, 18, 6]. Indeed, several secondary bifurcations appear, which include symmetry-breaking bifurcations and vortex-bearing states [10]. The bifurcation diagram plots the number of bosons in the BEC (1) as a function of the chemical potential . When , the non-linearity of the problem becomes irrelevant and the states bifurcate from the respective linear limit. Since an arbitrary potential can be approximated as a harmonic potential at the vicinity of a stable equilibrium point, when

we can decompose the linear eigenfunction

in Cartesian form as being proportional to

 |m,n⟩ := ϕm,n∼Hm(√Ωx)Hn(√Ωy)e−r22Ω,

where is the Hermite polynomial with being the associated quantum number of the harmonic oscillator. The eigenvalue corresponding to linear eigenfunction is . Notice that the characteristic eigenvalue parameter, i.e. the eigenvalue responsible for the bifurcation, is the chemical potential.

Let be a steady state solution to (4). In order to study its stability property we assume a perturbation ansatz around of the form

 ~Φ0(r,t)=e−iμt{ϕ0(r)+ϵ[a(r)eiωt+b(r)e−iω∗t]}, (5)

where is the complex eigenfrequency, is a “small” amplitude of the perturbation, and the star stands for complex conjugation. Plugging (5) into eq. (2), at order we obtain an eigenvalue problem that can be written in matrix form as:

 [A11A12−A∗12−A11]v=−ωv,v=[ab],

with eigenfrequencies

, and matrix blocks given by

 A1,1=−12Δ+2|ϕ0|2+W−μ,A1,2=|ϕ0|2.

The steady state is classified as stable in this Hamiltonian system if all the eigenfrequencies have vanishing imaginary part . On the other hand, when the solution becomes unstable, two types of instabilities can occur: i) exponential instabilities characterized by a pair of imaginary eigenfrequencies with zero real part , and ii) oscillatory instabilities characterized by a complex eigenfrequency quartet.

## 3 Numerical approximation of the problem

The Gross–Pitaevskii equation (2) describes a complex physical phenomenon with a rich bifurcating behavior, which is due to its nonlinearity. The bifurcation diagrams and the stability of the associated steady states can be studied numerically. The reconstruction of a bifurcation diagram entails computing multiple solution of a parametrized, non-linear problem, which can be computationally very demanding. For this reason we decided to take advantage of Reduced Order Methods (ROMs).

ROMs have been introduced for parametrized problems requiring real-time capabilities due to a many-query setting. The goal is to compute reliable results at a fraction of the cost of a conventional (full order) method. A practical way to realize this is to organize the computation in two steps:

• An offline phase: full order approximation solutions corresponding to selected representative parameters values/system configurations are computed and stored, together with other information concerning the parametrized problem. This is a computationally expensive step usually performed on high performance computing facilities.

• An online phase: the information obtained during the offline phase is used to compute the solution for a newly specified value of the parameters in a short amount of time (ideally in real time), even on a relatively low power device such as a laptop or a smartphone.

These split computational procedures are built in such a way that new parameter dependent quantities are easily and quickly computed online, while representative basis functions for selected parameter values and more demanding quantities are pre-computed offline. Among all the possible ROMs, we choose the Reduced Basis technique

The algorithm we use to rapidly plot the bifurcation diagram, presented in [20], is based on three building blocks :

1. A continuation technique: it permits to reconstruct properly the bifurcation diagram by following each branch and providing a proper initial guess for the non-linear iterations (next building block). We use a slight modification of the continuation method in [3], which consists in a for loop in the parameter domain where at each new cycle we check if the problem bifurcates. If it does bifurcate, we set that bifurcated solution as the initial guess for the next cycle.

2. Newton’s method: each point in the bifurcation diagram requires the solution of nonlinear eq. (4) for the corresponding value . To deal with the nonlinearity, we use the Newton-Kantorovich method [7].

3. A Galerkin Finite Element discretization: at each step of the Newton iteration, we have to approximate the solution of a new linear weak formulation. For the space discretization, we use the Galerkin-Finite Element method. This choice is made also in view of the numerical extension towards the model order reduction.

Algorithm 1 summarizes the steps for a generic nonlinear parametric problem that reads: given a scalar parameter , find , being a given functional space, such that

 G(X(μ);μ)=0. (6)

The next subsection is devoted to the description of building blocks 2 and 3.

### 3.1 Galerkin Finite Element method

As mentioned in the previous subsection, for the space discretization of the Gross–Pitaevskii equation we apply the Galerkin finite element method. Here, we introduce some standard notion for the discretization of generic problem (6) and then write the discrete counterpart of eq. (4).

In order to write the variational formulation the Galerkin Finite Element method builds upon, we introduce a family of finite dimensional spaces , such that . Let . Given , we seek that satisfies

 ⟨G(Xh(μ);μ),Yh⟩≐g(Xh(μ),Yh;μ)=0 ,∀ Yh∈Vh. (7)

We can not directly apply the finite element method to the equation (7) because of the non-linearity in . So, we first apply the Newton-Kantorovich method [7, 21], which reads as follows. Let be an initial guess. For every , we seek the variation such that

 dg[Xkh(μ)](δXh,Yh;μ)=g(Xkh(μ),Yh;μ) ,∀ Yh∈Vh, (8)

and then we update the solution for the successive step as , until we reach convergence for a prescribed tolerance .

From the algebraic point of view we denote with a basis for , such that the Newton’s method combined with the Galerkin finite element method and applied to our weak formulation reads: find such that

 J(→Xkh(μ);μ)δ→Xh=Gh(→Xkh(μ);μ) , (9)

where the Jacobian matrix in is defined as

 J(→Xkh(μ);μ))ij=dg[Xkh(μ)](Ej,Ei;μ),for alli,j=1,…,Nh. (10)

We can now present the weak formulation projected into the Finite Element space for problem (4). We recall that the solution to eq. (4) is a complex function. Let and be its real and imaginary part, respectively. The generic -th iteration of the Newton’s method (9) reads: seek , such that

 a(δXh,Yh)+b(δXh,Yh)−d(δXh,Yh;μ)+c(δXh,Xkh,Yh)= (11) a(Xkh,Yh)+b(Xkh,Yh)−d(Xkh,Yh;μ)+n(Yh)∀Yh∈Vh

where

 a(X,Y)=12∫D∇X⋅∇Ydx,b(X,Y)=12Ω∫D|r|2X⋅Ydx, d(X,Y;μ)=μ∫DX⋅Ydx,n(Y,Z)=∫D|Z|2Z⋅Ydx, c(X,Y,Z)=∫D[2(X⋅Z)Z+|Z|2X]⋅Ydx,

### 3.2 Reduced Basis method

In the previous section we derived a parametrized weak formulation that, once projected in a suitable Finite Element space, gives us a huge non-linear system which has to be solved for every parameter . Thus, the computational cost needed in order to study the bifurcating phenomena of the Gross-Pitaevskii equation has to be controlled, in such a way that we can explore rapidly and reliably all the branches. A possible strategy is the application of the reduced order methods (ROM), a collection of methodologies used to replace the original high dimension problem, typically called high fidelity approximation, with a reduced order one that is easy to manage.

One of these methodologies is the Reduced Basis method (RB), that roughly speaking consists in a projection of the high fidelity problem on a subspace of smaller dimension, constructed with some properly chosen basis functions.

As we have seen, the preliminary step is the projection of the weak formulation in a discrete space, which results in the Galerkin problem (7). Finding a numerical solution to this problem is very challenging because of the potential high number of degrees of freedom . Thus, we aim at constructing a discrete reduced order manifold, in which the parametrized solutions lie, that efficiently approximate the high fidelity one. This is the description of the first step, namely the off-line phase, in which we explore the parameter space in order to construct a basis for the low dimensional manifold, solving times the Galerkin high-fidelity problem varying .

On the other side, the second step, called on-line phase, is the efficient and reliable part where the solution is computed through the projection on the reduced manifold. In this way we do not have to project the equation in the huge Finite Element space for every , rather we can project on the low dimensional space spanned by the selected base. This complexity reduction is based on two main key points: the assumption that it holds the affine decomposition and the fact that . The latter assumption means that we can approximate the discrete manifold with a low number of basis functions.

To be precise we want to construct the reduced problem through the projection on a subspace spanned by a collection of the so called snapshots, i.e. the solution of the full order problem for different values of , obtained by Proper orthogonal decomposition (POD) or Greedy techniques [11, 19, 21].

As well as in the off-line phase (7), given , we seek that satisfies

 g(XN(μ),YN;μ)=0 ,∀ YN∈VN. (12)

At this point we have again to face up with the non-linearity, and so come back to the Newton-Kantorovich method obtaining: given an initial guess , for every we seek the variation such that

 dg[XkN(μ)](δXN,YN;μ)=g(XkN(μ),YN;μ) ,∀ YN∈VN, (13)

and then we update the solution as until convergence.

During the off-line phase, through the POD sampling, we obtain an orthonormal base given by , such that if we denote we can write every element as

 XN(μ)=N∑m=1X(m)N(μ)Σm , (14)

denoting with

the reduced solution vector.

Now we can choose the test element as for every , obtaining the algebraic system in given by

 g(N∑m=1X(m)N(μ)Σm,Σn;μ)=0 ,∀n=1,…,N . (15)

Moreover, we denote with

 (GN(→XN(μ);μ))n=g(N∑m=1X(m)N(μ)Σm,Σn;μ) ,

the residual reduced vector and with the transformation matrix whose elements are the nodal evaluation of the m-th basis function at the j-th node. Furthermore, with this new notation we observe that (15) corresponds to the solution of

 VTGN(V→XN(μ);μ)=0 . (16)

We can now recovery the algebraic low dimensional system by combining the Newton’s method and the Reduced Basis technique. Thus the problem reads as follows: find such that

 JN(→XkN(μ);μ)δ→XN=GN(→XkN(μ);μ) , (17)

where is the reduced Jacobian matrix defined as

 JN(→XkN(μ);μ)=VTJ(V→XkN(μ);μ)V . (18)

Furthermore, once finished the off-line phase, we have a base for the space , given by the span of the orthonormalized (through the Gram-Schmidt procedure) POD modes. In order to construct properly the bifurcation plots, the training set for the POD technique, with the values for which we solved the off-line problem, was setted as an ordered sampling of the interval .

Before ending this section we want to analyze in a more detailed way the system we obtained. Now we show how the Reduced Basis method reflects the projection properties of the Galerkin method also in the on-line phase. In fact, we can project the weak formulation (7) in the reduced manifold , obtaining the following: given and an initial guess for until convergence we seek such that

 a(δXN,YN)+b(δXN,YN)−d(δXN,YN;μ)+c(δXN,YN,XKN)= (19) a(XkN,YN)+b(XkN,YN)−d(XkN,YN;μ)+n(YN,XkN)∀YN∈VN

and then set . The matrix formulation follows directly and we obtain exactly the equation (17), where the reduced Jacobian is given by

 JN(→XkN(μ);μ)=AN+BN−μDN+CN (20)

where we have defined the reduced matrices as:

 AN=VTAhV ,BN=VTBhV , DN=VTDhV ,CN=N∑n=1X(n)NVTCh(Σn)V

in terms of the high fidelity ones coming from (17).

We have illustrated the on-line phase, that permits an efficient evaluation of the solution and possibly related outputs for every possible choice of a different parameter . The key point of this time saving is the so called affine decomposition. Indeed we want the computations to be independent form the, usually very high, number of degrees of freedom of the true discrete problem. In general the reduced matrices we have just presented are

-dependent and an affine-recovery technique called Empirical Interpolation Method (EIM) is needed

[5].

## 4 Results

All the results reported in this section have been obtained with FEniCS [1, 15, 4] for the high fidelity phase, and RBniCS [22] for the reduced order one.

### 4.1 Validation of the full order method

To validate the Full Order Method described in Sec. 3.1, we use a test proposed in [6]. We approximate the solution to eq. (4) in domain . Homogeneous Dirichlet boundary conditions are imposed on the entire boundary of . We recall that we set . For the space discretization, we use finite elements and a mesh with 6889 elements.

Fig. 1 displays the bifurcation diagram in the plane (left) and in the plane (right) obtained with the Full Order Method. These diagrams show the first three bifurcation points and the relative non-uniqueness of the solution with respect to the parameter . As is increased, the sequence of events observed in Fig. 1 is as follows. The ground state is the system simplest state. Its linear eigenfunction has corresponding eigenvalue . The ground state is generically stable, thus no further bifurcations occur from this state [12]. As expected, a unique solution branch departs from in Fig. 1. A representative density function for this branch is shown in Fig. 1(a). We see no further bifurcation for . The first interesting events in terms of bifurcation analysis occur for , with : two branches, associated to and , bifurcate from point in the and planes [17, 9]. Indeed, from point in Fig. 1 we observe the two expected branches. Representative density functions for these two branches are reported in Fig. 1(b) and 1(c). The next, more complicated, case of bifurcations emanates from point , with . In Fig. 1, we see that three branches depart from this point, associated to , , and . The corresponding representative densities are shown in Fig. 1(f), 1(e), and 1(d).

The density functions reported in Fig. 2, which are associated to and all the 6 solution branches in Fig. 1, show the richness of density patterns in order of decreasing In particular, we see the ground state in Fig. 1(a), the single charge vortex in Fig. 1(b), the 1-dark soliton stripe in Fig. 1(c), the dark soliton cross in Fig. 1(f), the ring dark soliton in Fig. 1(e), the 2-dark soliton stripe in Fig. 1(d). Notice that the 6 branches Fig. 1 are related to the first three eigenvalues. For example, the second bifurcation stems from a double eigenvalue and thus we have two branches, this phenomena is called multiple bifurcations. The stability property of these branches are different for each case, i.e the single charge vortex is always stable while the 1-dark soliton stripe is subject to multiple secondary bifurcations. These properties can be easily studied using standard techniques and are not among the purpose of this work.

To obtain the 6 branches of bifurcation diagrams in Fig. 1, we used a simple continuation method on parameter , combined with the aforementioned Finite Element discretization technique and the Newton’s method as non-linear solver. The overall simulation time required to complete the diagrams is roughly 95 minutes with continuation step . Consider that the mesh we used is quite coarse.

### 4.2 Reduced order approach

Next, we present the results obtained with our reduced order method as described in Sec. 3.2. Fig. 3 shows the first 6 bifurcation branches. We see that the curves obtained with the reduced order method matched very closely the curves computed with the full order method shown in Fig. 1. As further evidence of the accuracy of our ROM approach, we plot in Fig. 4 the difference between ROM and FOM solution for the solution and the density function in the and norms for each branch in Fig. 3.

As we said the result are exactly the same, but we gain a speedup of almost 1.1, so to construct this plots we took 5147 seconds. In fact as we can see in the next figures for the different branch, to obtain this computational advantage we paid a negligible worsening of the accuracy.

Now we can observe that we saved almost 600 seconds, and still our computations are not independent from the number of degrees of freedom of the Finite Element Method. To improve a lot the results in terms of speed-up we need to implement some affine recovery technique, such as EIM or DEIM.

So we can present again the bifurcation diagrams and the representative solutions of each branch.

We investigated also the 2 parameters case, considering the trap strength as a parameter varying in . In this case so we want to detect the evolution of the first bifurcations for different trapping configuration.

### 4.3 Hyper-reduction techniques

We used EIM/DEIM for the construction of one of the branch of the second bifurcation and we obtained results comparable with the previous ones as we can see in the next figures.

Moreover, the error plots the results are exactly the same in terms of accuracy with respect to the affine recovery techniques.

In both this cases we obtained a good performance in terms of reduced error plot, but also we improved a lot the computational savings. In fact here the on-line phase is independent from the number of degrees of freedom of the finite element method, so the RB method is really efficient.

For the EIM we took just 55 seconds to construct the reduced basis bifurcation diagrams with respect to the 246 seconds for the high fidelity one, so obtaining a speed-up of almost 5. Moreover we obtained a surprisingly better result using DEIM. In fact it tooks just 7 seconds for the complete construction of the reduced basis bifurcation diagrams, with respect to the 246 seconds for the high fidelity one, so obtaining a speed-up of almost 32.

The hyper-reduction techniques are crucial in the two parameters test case. In fact, in order to reconstruct the diagram in Fig. 7 we need 14708 seconds while applying the DEIM technique we just need 366 seconds, with a speed-up of more then 40 times.

## Acknowledgments

This work was supported by European Union Funding for Research and Innovation through the European Research Council (project H2020 ERC CoG 2015 AROMA-CFD project 681447, P.I. Prof. G. Rozza). This work was also partially supported by NSF through grant DMS-1620384 (A. Quaini). We acknowledge fruitful conversations with Dr. Sean Garner.

## References

• [1] FEniCS Project.
• [2] G L Alfimov and D A Zezyulin. Nonlinear modes for the gross–pitaevskii equation—a demonstrative computation approach. Nonlinearity, 20(9):2075, 2007.
• [3] E. Allgower and K. Georg. Introduction to Numerical Continuation Methods. Society for Industrial and Applied Mathematics, 2003.
• [4] Martin S. Alnæs, Jan Blechta, Johan Hake, August Johansson, Benjamin Kehlet, Anders Logg, Chris Richardson, Johannes Ring, Marie E. Rognes, and Garth N. Wells. The fenics project version 1.5. Archive of Numerical Software, 3(100), 2015.
• [5] M. Barrault, N. C. Nguyen, Y. Maday, and A. T. Patera.

An “empirical interpolation” method: Application to efficient reduced-basis discretization of partial differential equations.

C. R. Acad. Sci. Paris, Série I., 339:667–672, 2004.
• [6] E.G. Charalampidis, P.G. Kevrekidis, and P.E. Farrell. Computing stationary solutions of the two-dimensional gross–pitaevskii equation with deflated continuation. Communications in Nonlinear Science and Numerical Simulation, 54:482 – 499, 2018.
• [7] P. G. Ciarlet. Linear and Nonlinear Functional Analysis with Applications:. Other Titles in Applied Mathematics. Society for Industrial and Applied Mathematics, 2013.
• [8] M P Coles, D E Pelinovsky, and P G Kevrekidis. Excited states in the large density limit: a variational approach. Nonlinearity, 23(8):1753, 2010.
• [9] Andres Contreras and Carlos Garcia-Azpeitia. Global bifurcation of vortex and dipole solutions in bose-einstein condensates. Comptes Rendus Mathematique, 354(3):265 – 269, 2016.
• [10] C. García-Azpeitia and D. E. Pelinovsky. Bifurcations of multi-vortex configurations in rotating bose–einstein condensates. Milan Journal of Mathematics, 85(2):331–367, Dec 2017.
• [11] J. Hesthaven, G. Rozza, and B. Stamm. Certified Reduced Basis Methods for Parametrized Partial Differential Equations. Springer Briefs in Mathematics, 2015.
• [12] P. Kevrekidis, D. Frantzeskakis, and R. Carretero-González. The Defocusing Nonlinear Schrödinger Equation. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2015.
• [13] P G Kevrekidis, V V Konotop, A Rodrigues, and D J Frantzeskakis. Dynamic generation of matter solitons from linear states via time-dependent scattering lengths. Journal of Physics B: Atomic, Molecular and Optical Physics, 38(8):1173, 2005.
• [14] Yuri S. Kivshar, Tristram J. Alexander, and Sergey K. Turitsyn. Nonlinear modes of a macroscopic quantum oscillator. Physics Letters A, 278(4):225 – 230, 2001.
• [15] Anders Logg, Kent-Andre Mardal, Garth N. Wells, et al. Automated Solution of Differential Equations by the Finite Element Method. Springer, 2012.
• [16] S. Middelkamp, P. G. Kevrekidis, D. J. Frantzeskakis, R. Carretero-González, and P. Schmelcher. Bifurcations, stability, and dynamics of multiple matter-wave vortex states. Phys. Rev. A, 82:013646, Jul 2010.
• [17] S. Middelkamp, P. G. Kevrekidis, D. J. Frantzeskakis, R. Carretero-González, and P. Schmelcher. Bifurcations, stability, and dynamics of multiple matter-wave vortex states. Phys. Rev. A, 82:013646, Jul 2010.
• [18] S. Middelkamp, P.G. Kevrekidis, D.J. Frantzeskakis, R. Carretero-González, and P. Schmelcher. Emergence and stability of vortex clusters in bose-einstein condensates: A bifurcation approach near the linear limit. Physica D: Nonlinear Phenomena, 240(18):1449 – 1459, 2011.
• [19] A.T. Patera and G. Rozza.

Reduced basis approximation and A posteriori error estimation for Parametrized Partial Differential Equation

.
MIT Pappalardo Monographs in Mechanical Engineering, Copyright MIT (2007-2010).
• [20] F. Pichi and G. Rozza. Reduced basis approaches for parametrized bifurcation problems held by nonlinear Von Kármán equations. Journal of Scientific Computing, 2019.
• [21] A. Quarteroni, A. Manzoni, and F. Negri. Reduced Basis Methods for Partial Differential Equations: An Introduction. UNITEXT. Springer International Publishing, 2015.
• [22] RBniCS.