Mathematics Subject Classification (2020): Primary 65M12; Secondary 65M15, 65L20, 41A30, 65D15

Keywords: parametric evolution equations in Hilbert spaces; nonlinear manifold approximation; residual-based time stepping; convergence analysis; Gaussian manifold

1 Introduction↩︎

The numerical solution of time-dependent partial differential equations in high dimensions remains a central challenge in scientific computing. In particular, evolution equations of Schrödinger type arise ubiquitously in quantum dynamics, molecular physics, and wave propagation, where the solution typically exhibits oscillatory behavior, strong localization, and multi-scale structure. Standard grid-based discretizations quickly become infeasible as the spatial dimension increases, whereas reduced-order approaches must carefully balance accuracy, stability, and computational cost.

We consider time-dependent evolution equations posed in the Hilbert space \(\mathcal{H}= L^2(\Omega,\mathbb{C})\) with a general right-hand side operator \(f:[0,T]\times \Omega \times \mathcal{H}\to \mathbb{C}\) of the form \[\label{eq:PDE95original} \left\{ \begin{array}{ll} \partial_t \psi(t,x) = f(t,x,\psi), & (t,x)\in [0,T]\times \Omega, \\[1ex] \psi(0,x) = \psi_0(x), & x\in \Omega . \end{array} \right.\tag{1}\] The initial datum \(\psi_0\in \mathcal{H}\) may, for instance, be normalized with \(\|\psi_0\|=1\). The right-hand side \(f(t,\cdot,\psi)\in \mathcal{H}\) may represent a time-dependent Schrödinger operator, \[f(t,x,\psi) = \frac{1}{{\rm i}}\left(-\frac{1}{2} \Delta_x \psi(t,x) + V(t,x)\psi(t,x)\right), \quad (t,x)\in [0,T]\times \Omega,\] or, more generally, a linear partial differential operator of the form \[f(t,x,\psi) = \sum_{\alpha} c_{\alpha}(t,x)\, \frac{\partial^{|\alpha|}}{\partial x^\alpha}\psi(t,x), \quad (t,x)\in [0,T]\times \Omega.\] Nonlinear evolution operators are likewise admissible, provided that the resulting initial value problem is well posed in \(\mathcal{H}\), and in particular allows for some type of (one-sided) Lipschitz control.

A powerful class of methods approximates the solution \(\psi(t,\cdot)\) within a low-dimensional nonlinear manifold \(\mathcal{M}\) embedded in \(\mathcal{H}\). Rather than resolving the full spatial domain, one represents the solution in the form \[\psi(t,\cdot) \approx u(\theta(t),\cdot),\] where the parameter vector \(\theta(t)\in\Theta\) evolves in time, for a complex parameter space \(\Theta\). Such parametrized ansatz spaces include Gaussian wave packets [1], [2], Hagedorn bases [3], [4], neural networks [5]–[8] also with Gaussian activation functions [9]–[11], and more general manifolds. These representations can capture localization, phase information, and deformation using comparatively few degrees of freedom.

A central mathematical framework underlying many parametrized evolution schemes is the Dirac–Frenkel variational principle [12], which projects the infinite-dimensional dynamics onto the tangent space of the approximation manifold. This yields a system of ordinary differential equations for the parameters \(\theta(t)\) obtained through a linear least-squares condition. While the continuous-time formulation enjoys attractive geometric properties, including conservation of norm and energy in Hamiltonian systems [13], its numerical implementation suffers from the so-called matrix singularity problem [14] or tangent space collapse [15]. An alternative viewpoint reverses the order of discretization and parametrization. One first discretizes the original evolution equation in time and subsequently parametrizes within the manifold each time step by minimizing a residual in the Hilbert space norm. This “discretize-then-parametrize” philosophy leads naturally to a sequence of nonlinear optimization problems for the parameters [16], [17]. Both approaches are actively used in practice and give rise to different numerical schemes with distinct stability and convergence properties. Despite their widespread use, a unified mathematical understanding of these parametric time-stepping strategies remains incomplete, and error bounds that account simultaneously for time discretization and manifold approximation seem to be rare. Notable exceptions are [15], which motivated the present work, the analysis of regularized schemes [18], [19], and [10] for a continuous-time analysis of dynamical sampling [10].

Here, we aim at a unified framework for the analysis of time-discrete parametrized approximations of evolution equations in Hilbert spaces. We consider nonlinear approximation manifolds of the form \[\mathcal{M} = \{ u(\theta,\cdot)\in\mathcal{H}\;: \;\theta \in \Theta \subset \mathbb{C}^p, p\geq 1 \},\] and formulate numerical schemes through the minimization of residual functions at each time step. Two natural residual structures arise: (i) discretization–then–parametrization residuals, obtained by first applying a time discretization scheme within the family of \(\zeta\)-methods and subsequently optimizing the parameters on the manifold; (ii) parametrization–then–discretization residuals, derived from time discretizations of the Dirac–Frenkel variational principle. Both approaches lead to iterative schemes of the form \[\begin{align} &\theta_{k+1} = \arg\min_{\theta\in\Theta} \| r_k(\theta,\cdot) \|_{\mathcal{H}},\\ &u(\theta_{k+1},\cdot) \approx \psi(t_k,\cdot),\quad k\ge 0, \end{align}\] where \(r_k(\theta,\cdot)\) denotes the corresponding residual at the \(k\)th time step.

Our first main contribution is a unified error analysis for both residual formulations without assumptions on the underlying approximation manifold, see 1 and 2. We show that the global approximation error decomposes into a residual contribution arising from the manifold optimization and a classical time discretization contribution governed by the regularity of the exact solution and the choice of \(\zeta\)-scheme. This yields explicit convergence rates of order \(h\) or \(h^2\) under natural Lipschitz, one-sided Lipschitz, and dissipativity assumptions. For the parametrization–then–discretization approach, additional conditioning terms involving the singular values of the parametrizing gradient appear naturally and clarify stability requirements. A second contribution concerns the explicit evaluation of residual norms and their gradients for Gaussian approximation manifolds. When the evolution operator acts polynomially on Gaussian functions, all inner products required for residual minimization admit closed-form expressions in terms of Gaussian moments. We present compact formulas for these integrals and their parameter derivatives, which pave the way to efficient computation even in higher dimensions. When closed forms are unavailable, we briefly discuss quadrature-based approximations and their accuracy. This analytic tractability highlights an important advantage of Gaussian manifolds: spatial discretization can be avoided to a large extent, and the computational effort is concentrated on low-dimensional optimization problems.

To illustrate the theoretical results, we present numerical experiments for one-dimensional Schrödinger equations with harmonic, double-well, and hyperbolic cosine potentials. These examples probe regimes of exact representability, dynamical interference, and oscillatory wave packet motion. The experiments confirm the predicted convergence orders and demonstrate how residual minimization accuracy directly influences accuracy and conservation behavior. Although our numerical investigations focus on Schrödinger dynamics, the framework is not restricted to conservative evolution equations. To illustrate this broader applicability, Appendix 12 briefly discusses a Fokker–Planck example.

More generally, the framework developed here connects variational approximation theory, time discretization of evolution equations, and nonlinear optimization within a single analytical structure. While Gaussian manifolds and time-dependent Schrödinger equations serve as a concrete and computationally attractive example, the error analysis applies to general nonlinear parametrizations and evolution equations.

Structure of the paper↩︎

2 introduces the abstract parametrized approximation framework and the residual minimization formulation. 3 derives explicit residuals corresponding to discretize-then and parametrize-then strategies within the family of \(\zeta\)-methods. 4 develops the unified error analysis. 5 presents closed-form Gaussian integral formulas, quadrature approximations for residual evaluation, and details of the initialization scheme. 6 illustrates the dynamics of the chosen examples. 7 reports numerical experiments illustrating the theoretical results. Conclusions are given in 8.

2 Parametrized dynamical approximation↩︎

We work in an ambient Hilbert space \(\mathcal{H}= L^2(\Omega,\mathbb{C})\) with \(\Omega\subseteq\mathbb{R}^d\) and denote the inner product and induced norm by \(\langle\cdot,\cdot\rangle_\mathcal{H}\) and \(\Vert\cdot\Vert_\mathcal{H}\), respectively. For solving the evolutionary equation 1 , we introduce a time grid \(t_k = k h\) of uniform step-size \(h>0\) and aim at an optimal parameter sequence \((\theta_k)_{k\ge0}\subset\Theta\) from the parameter set \(\Theta\subset\mathbb{C}^p\) with \(p>1\). We aim at the approximation of the PDE solution, \(\psi(t_k,x)\approx u(\theta_k,x)\) for all time steps \(k\ge0\). Given a residual function \[\begin{align} r_k: \Theta\times\Omega \to \mathbb{C}, \quad (\theta,x) \mapsto r_k(\theta,x), \end{align}\] parametrized by a ‘known’ parameter \(\theta_k\) at time step \(t_k\), we define the cost function \[F_k: \Theta\subset\mathbb{C}^p \to \mathbb{R}, \quad \theta \mapsto \Vert r_k(\theta,\cdot) \Vert^2_{\mathcal{H}}.\] This is the function to be minimized in the parameter space \(\Theta\subset\mathbb{C}^p\) to determine the next point in the discrete trajectory of approximate solutions. Namely, given an initial condition \(\psi_0\in\mathcal{H}\), we determine the approximation \(u(\theta_0,\cdot)\in\mathcal{M}\) and then compute the subsequent values \(u(\theta_{k+1},\cdot)\in\mathcal{M}\) iteratively via the optimization sequence \[\theta_{k+1} = \arg \min_{\theta\in\Theta} F_k(\theta),\quad k=0,1,\ldots\] We will consider and analyze two specific choices for the residual function \(r_k(\theta,\cdot)\) later in 3.

2.1 Abstract discretization process↩︎

Given a choice of step-size \(h>0\), the required number of steps \(N\) in the interval \([0,T]\) is given by \(N=\lceil T/h \rceil\). For the time grid \(t_0,\ldots,t_N\), we perform the following process to determine the discrete solutions \(u(\theta_0,\cdot),\ldots,u(\theta_N,\cdot)\):

• Given the initial condition \(\psi_0\in\mathcal{H}\), obtain \(u_0\in\mathcal{M}\) via \[\theta_{0} = \arg\min_{\theta\in\Theta} \Vert \psi_0(\cdot) - u(\theta,\cdot) \Vert_\mathcal{H}\]

• Construct the cost function \(F_k(\theta)=\Vert r_k(\theta,\cdot)\Vert_\mathcal{H}^2\) for \(k=0\).

• Minimize \(F_k\) and assign the optimum to \(\theta_{k+1}\) for \(k=0\).

• Repeat steps (ii) and (iii) for \(k=1,2,\ldots\) to obtain \(\theta_2,\theta_3,\ldots\)

• Compute \(u(\theta_k,\cdot)\) at the values \(\theta_k\) determined in steps (i)–(iv) to provide the discrete approximations \(\psi(t_k,\cdot)\approx u(\theta_k,\cdot)\) for all \(k=0,1,\ldots,N\)

See also 1 for a more algorithmic formulation of the process. Our main focus here will be on the comparative analysis of several natural choices of residual optimization, which result from different orderings of time discretization and parameter optimization, as well as a different choice of the time-integration within the class of \(\zeta\)-methods.

2.2 Initialization↩︎

The first step in the entire process is determining the appropriate parameter \(\theta_0\in\Theta\) that approximates the initial condition \(\psi_0\in\mathcal{H}\). To do this, we define the cost function \[\label{eq:cost95Fn0} \begin{align} F(\theta) &= \Vert \psi_0 - u(\theta,\cdot) \Vert^2_\mathcal{H}, \\ &= \langle \psi_0,\psi_0 \rangle_\mathcal{H}- 2 \Re \langle\psi_0, u(\theta,\cdot)\rangle_\mathcal{H}+ \langle u(\theta,\cdot),u(\theta,\cdot)\rangle_\mathcal{H}. \end{align}\tag{2}\] For example, we might consider a univariate mixture \[u(\theta,x) = \sum_{l=1}^L b(\theta^l,x)\] with holomorphically parametrized complex Gaussians \[\label{eq:gaussian95basis} b(\theta^l,x) = \exp(-a^l x^2 + \kappa^l x + \gamma^l),\quad x\in\mathbb{R},\tag{3}\]

with \(a\in \mathbb{H}= \{z\in\mathbb{C}: \Re(z) >0\}\), \(\kappa\in\mathbb{C}\), \(\gamma\in\mathbb{C}\) and total parameter \(\theta=(a,\kappa,\gamma)\), see also 5. Then, the gradient of the cost function satisfies \[\begin{align} \nabla_{\theta} F(\theta) = - \langle \psi_0,\nabla_\theta u(\theta,\cdot) \rangle_\mathcal{H}+ \langle u(\theta,\cdot),\nabla_\theta u(\theta,\cdot) \rangle_\mathcal{H}. \end{align}\] In this specific situation, we employ a heuristic five-step strategy to determine a minimizer \(\theta_0 = \arg\min_{\theta\in\Theta}F(\theta)\) with \(\theta_0=(a_0^l,\kappa_0^l,\gamma_0^l)_{l=1,\ldots,L}\), see 5.3 for further details.

2.3 Optimization Process↩︎

At each time step, an optimization problem is solved. The parameter \(\theta_{k+1}\) is defined as a minimizer of \(F_k\), namely, \[\theta_{k+1} = \arg\min_{\theta\in\Theta} F_k(\theta) ,\] where we recall \(F_k(\theta)=\Vert r_k(\theta,\cdot) \Vert^2_{\mathcal{H}}\). Since a global minimizer may not exist or may be computationally inaccessible, in particular with a constraint like that for the width matrix \(A\in\mathbb{H}^d\) of a Gaussian basis function, we only expect to obtain a parameter value \(\theta_{k+1}\) for which the residual is sufficiently small in a small neighborhood of the starting point \(\theta_k\) for every step \(k\). We emphasize that local minimizers with small enough residual will be our target here, since the overall approximation scheme has inherent time discretization, initialization, and quadrature errors. For the numerical experiments presented in 6, we perform unconstrained optimization in \(\mathbb{C}^p\) using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. This algorithm searches for the local minimum until one of the following stopping criteria is met: \(\Vert \nabla F_k(\theta) \Vert_{\infty} < \varepsilon_1\), \(|F_k(\theta^{(n)})-F_k(\theta^{(n-1)})|<\varepsilon_2\) or \(\Vert \theta^{(n)}-\theta^{(n-1)}\Vert_2 < \varepsilon_3\) for \(n=1,2,\ldots\), given predefined accuracies \(\varepsilon_1,\varepsilon_2,\varepsilon_3>0\). If none of these criteria is met, we terminate the computations if a certain maximum number of steps \(n_{\max}\) is reached. This method is available in various optimization packages, and the solver’s efficiency is significantly improved (especially in higher dimensions) if the explicit form of the gradient \(\nabla_\theta F_k(\theta)\) is provided. We discuss such an explicit form in 3.3.

The numerical experiments shown in [sec:Num_Experiments,sec:Num_Experiments_II] were carried out using unconstrained optimization. However, since these examples deal with Gaussian bases, one may impose the positivity of the real part of the width parameter \(a\) of each Gaussian, i.e., we choose a small \(\epsilon\) and impose the linear constraint \(a>\epsilon\). In our experiments, it was not necessary.

3 Residual and Cost Function↩︎

We consider a numerical approximation to the PDE 1 via a \(\theta\)-parametrization and an iterative process in the parameter space \(\Theta\subset \mathbb{C}^p\), where \(p\) is the complex dimension. The set of time-discrete approximations \(u(\theta_0,\cdot),u(\theta_1,\cdot),\ldots\) are constructed via the iterative process \[\label{eq:iteration} \theta_{k+1} = \arg \min_{\theta\in\Theta} \Vert r_k(\theta,\cdot) \Vert_{\mathcal{H}}, \quad k=0,1,\ldots\tag{4}\] for various forms of the residual function \(r_k(\theta,\cdot)\), whose explicit form will depend on the chosen numerical schemes.

3.1 Discretization then parametrization↩︎

The discretize-then approach falls into the class of Rothe methods for evolution equations, see for example the classic monograph [20]. In the context of parametric approximation for Schrödinger evolution, it has been proposed in [16], [17], [21]; see also [19] for more general stiff evolution equations. Here, following [15], we employ the usual \(\zeta\)-scheme for time-discretization, which generalizes the Euler methods and the trapezoidal rule. Namely, given a continuous function \(\varphi:[0,T]\to \mathbb{R}\), the integral over small intervals \([t_a,t_b]\) with length \(h=t_b-t_a\) are approximated via the weighted rule \(\int_{t_a}^{t_b} \varphi(t) {\rm d}t \approx h ( \zeta \varphi(t_a) + \hat{\zeta} \varphi(t_b))\) with \(\zeta\in [0,1]\) and \(\hat{\zeta}=1-\zeta\). The construction of the corresponding residual proceeds in two steps.
Step (i). Integrate both sides of the PDE 1 with respect to \(t\) on small intervals of the form \([t_{k+1},t_k]\) with length \(h=t_{k+1}-t_k\). Namely, \[\begin{align} \psi(t_{k+1}) - \psi(t_k) &= \int_{t_k}^{t_{k+1}} f(t,\cdot,\psi) {\rm d}t, \\ \psi_{k+1} - \psi_k &\approx h \left[ \zeta f(t_k,\cdot,\psi_k) + \hat{\zeta} f(t_{k+1},\cdot,\psi_{k+1}) \right] . \end{align}\] This leads to an abstract optimization problem in the Hilbert space \(\mathcal{H}\), \[\psi_{k+1} = \arg\min_{\psi\in\mathcal{H}} \Vert \psi - \psi_k - h [\zeta f(t_k,\cdot,\psi_k) + \hat{\zeta} f(t_{k+1},\cdot,\psi) ] \Vert_{\mathcal{H}} .\] Step (ii). Consider the \(\theta\)-parametrization \(\psi(t,\cdot)\approx u(\theta(t),\cdot)\) for all \(t\in [0,T]\) and replace these values in the expression shown above. This leads to an optimization problem in the parameter space \(\Theta\subset\mathbb{C}^p\), and we define the Discretization then Parametrization approach as the iterative method \[\theta_{k+1} = \arg \min_{\theta\in\Theta} \Vert r_k(\theta,\cdot) \Vert_{\mathcal{H}}, \quad k=0,1,\ldots\] with the residual function \[\label{eq:residual95Disc-Opt} r_k(\theta,\cdot) = u(\theta,\cdot) - u(\theta_k,\cdot) - h \left[ \zeta f(t_k,\cdot,u(\theta_k,\cdot)) + \hat{\zeta}f(t,\cdot,u(\theta,\cdot)) \right],\tag{5}\] see also [15]. Given the parameter \(\theta_k\in\Theta\) at time \(t_k\), we have in mind that an optimized parameter \(\theta_{k+1}\) provides an approximation \(u(\theta_{k+1},\cdot)\approx\psi(t_{k+1},\cdot)\). The value \(\zeta=1\) (resp. \(\zeta=0\)) corresponds to the forward (resp. backward) Euler method.

3.2 Parametrization then discretization↩︎

The second approach parametrizes first and then discretizes in time. The formal construction follows three steps. First, we consider the general Dirac-Frenkel approximation at a point \(u\in\mathcal{M}\). Second, we parametrize \(u\) with \(\theta\in\Theta\). Finally, we discretize in time considering a \(\zeta\)-scheme.
Step (i). Let \(\mathcal{M}\) be the manifold approximation to the space of solutions of the PDE 1 . Then, following the Dirac-Frenkel principle, we project \(f\) to the tangent space at \(u\in\mathcal{M}\), \[\partial_t u = P_{T_u \mathcal{M}} f(t,\cdot,u),\] where the projection is determined as the solution of the linear least squares problem \[P_{T_u \mathcal{M}} f(t,\cdot,u) = \arg\min_{w\in T_u\mathcal{M}} \Vert f(t,\cdot,u) - w \Vert_{\mathcal{H}} .\] Step (ii). Consider the parametrized manifold \(\mathcal{M}=\{ u(\theta,\cdot) \;: \;\theta\in\Theta \}\). This implies that \(\partial_t u = \nabla_{\theta} u \cdot \dot{\theta}\) and that the tangent space at \(u\in\mathcal{M}\) is of the form \(T_u \mathcal{M}= \{ \nabla_{\theta} u \cdot \eta \;: \;\eta\in\mathbb{C}^p \}\), leading to \[\label{eq:DF} \dot{\theta} = \arg \min_{\eta\in\mathbb{C}^p} \Vert \tilde{r}(\theta,\eta,\cdot) \Vert_{\mathcal{H}} ,\tag{6}\] where \(\tilde{r}(\theta,\eta,\cdot)= \nabla_{\theta} u \cdot \eta - f(t,\cdot,u)\) is the Dirac–Frenkel residual function.
Step (iii). Consider time discretization via the approximation \(\dot{\theta}|_{t_k} \approx (\theta_{k+1}-\theta_k)/h\). Thus, the tangent space at \(u_k\) is constrained to \(\eta=(\theta-\theta_k)/h\), i.e., \(T_{u_k} \mathcal{M}= \{ \nabla_{\theta} u_k \cdot (\theta-\theta_k)/h \;: \;\theta\in\Theta \}\), leading to \[\theta_{k+1} = \arg\min_{\theta \in \Theta} \Vert \nabla_{\theta_k} u_k \cdot (\theta-\theta_k) - h f(t,\cdot,u_k) \Vert_{\mathcal{H}} .\] Similarly, the discretization \(\dot{\theta}|_{t_{k}} \approx (\theta_k-\theta_{k-1})/h\) (after a relabelling \(k\leftarrow k+1\)) leads to a search inside of the moving tangential space \(V_u = \{ \nabla_{\theta} u \cdot (\theta-\theta_k)/h \;: \;\theta\in\Theta \}\), where \(u=u(\theta,\cdot)\). Thus, \[\theta_{k+1} = \arg\min_{\theta \in \Theta} \Vert \nabla_{\theta} u \cdot (\theta-\theta_k) - h f(t,\cdot,u) \Vert_{\mathcal{H}} .\] The first expression corresponds to the forward Euler method, whereas the second corresponds to the backward Euler method. Consider an extension to \(\zeta\in [0,1]\) such that \(\zeta=1\) (\(\zeta=0\)) will reduce to forward (backward) Euler. This condition is satisfied by the residual function \[\label{eq:residual95Opt-Disc} r_k(\theta,\cdot) = \left[ \zeta \nabla_{\theta_{k}} u(\theta_k,\cdot) + \hat{\zeta} \nabla_{\theta} u(\theta,\cdot) \right] \cdot (\theta-\theta_k) - h \left[ \zeta f(t,\cdot,u(\theta_k,\cdot)) + \hat{\zeta} f(t,\cdot,u(\theta,\cdot)) \right] .\tag{7}\]

Note that in comparison to the Dirac–Frenkel residual \(\tilde{r}(\theta,\dot{\theta},\cdot)\), there is a multiplicative factor \(h\), namely, \(r_k(\theta_{k+1},\cdot) \approx h \tilde{r}(\theta|_{t_k},\dot{\theta} |_{t_k}, \cdot)\). More precisely, \[\lim_{h\to 0} r_k(\theta_{k+1},\cdot)/h = \tilde{r}(\theta|_{t_k},\dot{\theta} |_{t_k}, \cdot) .\]

3.2.1 Euler methods↩︎

For the two extreme cases \(\zeta\in\{0,1\}\), the residual 7 has been considered before. For \(\zeta=1\), it takes the form \[r_k(\theta_{k+1},\cdot) = \nabla_\theta u(\theta_k,\cdot) (\theta_{k+1} - \theta_k) - h f(t_k,\cdot,u(\theta_k,\cdot)) .\] The corresponding one-step method \(\theta_{k+1} = \arg\min_{\theta\in\Theta}\|r_k(\theta,\cdot)\|_\mathcal{H}\) is the forward Euler method for the Dirac–Frenkel principle 6 , that is, \(\theta_{k+1} = \theta_k + h\dot{\theta}_k\) with \[\dot{\theta}_k = \arg \min_{\eta\in\mathbb{C}^p} \Vert \nabla_\theta u(\theta_k,\cdot) \eta - f(t_k,\cdot,u(\theta_k,\cdot)) \Vert_\mathcal{H}.\] For each time step, it requires to solve a linear least squares problem. On the other hand, for \(\zeta=0\), we obtain \[r_k(\theta_{k+1},\cdot) = \nabla_\theta u(\theta_{k+1},\cdot) (\theta_{k+1} - \theta_k) - h f(t_{k+1},\cdot,u(\theta_{k+1},\cdot)),\] see also [15] and for a regularized implicit Euler discretization [18]. Here, in the fully implicit approach, a nonlinear optimization is performed for each time-step. Note, that the implicit Euler method applied to the normal equations of the Dirac–Frenkel principle 6 , \[\langle \nabla_\theta u(\theta,\cdot),\nabla_\theta u(\theta,\cdot)\rangle_\mathcal{H}\;\dot{\theta} = \langle \nabla_\theta u(\theta,\cdot),f(t,\cdot,u(\theta,\cdot))\rangle_\mathcal{H},\] results in the related, but different nonlinear system \[\theta_{k+1} = \theta_k + h \langle \nabla_\theta u(\theta_{k+1},\cdot),\nabla_\theta u(\theta_{k+1},\cdot)\rangle_\mathcal{H}^{-1}\langle \nabla_\theta u(\theta_{k+1},\cdot),f(t_{k+1},\cdot,u(\theta_{k+1},\cdot))\rangle_\mathcal{H}.\]

3.3 Explicit form of the cost function↩︎

In this subsection and for the numerical experiments in [sec:Num_Experiments,sec:Num_Experiments_II], we consider a specific type of parametric approximation \(\psi(t,x)\approx u(\theta(t),x)\), namely, \[\label{eq:sum} u(\theta(t),x) = \sum_{l=1}^L b(\theta^l(t),x), \qquad x\in\Omega,\tag{8}\] where \(\theta(t)=(\theta^1(t),\ldots,\theta^L(t))\) with \(\theta^l(t)\in\Theta\) and \(b(\theta^l(t),\cdot)\in\mathcal{H}= L^2(\Omega,\mathbb{C})\) are the individual basis functions. We assume that the function \(f\) defining the PDE 1 acts on such an ansatz \(u\) as \[\label{eq:f95condition95L} f(t,x,u) = \sum_{l=1}^L \phi(\theta^l(t),t,x) b(\theta^l(t),x),\tag{9}\] with a function \(\phi\) such that \(\phi(\theta^l(t),t,\cdot)b(\theta^l(t),\cdot) \in \mathcal{H}\), for all \(l=1,\ldots,L\) and for all \(t\in [0,T]\). The following remark illustrates the above assumptions and demonstrates their natural occurrence in a quantum application.

3.3.1 Discretize-then↩︎

Assuming the basis form 8 and the \(f\)-condition 9 , the residual function 5 for the discretize-then approach is the finite sum \[\begin{align} \label{eq:residual95sum} r_k(\theta,\cdot) = \sum_{l=1}^L r^l_k(\cdot) \end{align}\tag{10}\] with summands \(r^l_k(\cdot) = Q_1(\theta^l_k,\theta^l,\cdot) b(\theta^l_k,\cdot) + Q_2(\theta^l_k,\theta^l,\cdot) b(\theta^l,\cdot)\), that depend on the basis functions \(b(\theta^l_k,\cdot)\) and \(b(\theta^l,\cdot)\) as well as the functions \[\begin{align} Q_1(\theta^l_k,\theta^l,\cdot) &= -1 - h \zeta \phi(\theta^l_k,\cdot) , \\ Q_2(\theta^l_k,\theta^l,\cdot) &= 1 - h \hat{\zeta} \phi(\theta^l,\cdot) . \end{align}\]

3.3.2 Parametrize-then↩︎

To express the parametrize-then residual 7 in an analogous form, we need an additional assumption. We have to require the existence of a vector \(v:\Omega\to\mathbb{R}^p\) such that \[\label{eq:grad95b} \nabla_{\theta} b(\theta,x) = \vec{v}(x) b(\theta,x),\tag{11}\] where the vector \(\vec{v}(x) = \vec{v}(\theta,x)\) might depend on \(\theta\) as well. For instance, \(\vec{v}(x) = (-x^2,x,1)\) if the basis function \(b(\theta,x)\) is the one-dimensional Gaussian 3 . In higher dimensions \(d\ge 1\), \(\vec{v}(x)\cdot\theta\) would correspond to the argument of the exponential, i.e., \(\vec{v}(x)\cdot\theta=-x^\top A x + \kappa^\top\cdot x + \gamma\). Then, we also have 10 with summands defined by \[\begin{align} Q_1(\theta^l_k,\theta^l,\cdot) &= \zeta \left[ \vec{v}(x)\cdot (\theta^l-\theta^l_k) - h \phi(\theta^l_k,\cdot) \right], \\ Q_2(\theta^l_k,\theta^l,\cdot) &= \hat{\zeta} \left[ \vec{v}(x)\cdot (\theta^l-\theta^l_k) - h \phi(\theta^l,\cdot) \right]. \end{align}\]

3.3.3 Unified form↩︎

Hence, for a parametrization \(u(\theta,\cdot) = \sum_{l=1}^L b(\theta^l,\cdot)\) with properties 9 and potentially 11 , the cost function for both residual approaches is explicitly given by \[\begin{align} F_k(\theta) = \Vert r_k(\theta,\cdot)\Vert_\mathcal{H}^2 = \sum_{l',l=1}^L \langle r_k^{l'}(\cdot) , r_k^l(\cdot) \rangle_\mathcal{H}, \end{align}\] with \[\begin{align} \langle r_k^{l'}(\cdot) , r_k^l(\cdot) \rangle_\mathcal{H}= \sum_{\alpha,\beta=1}^2 \langle Q_\alpha(\theta_k^{l'},\theta^{l'},\cdot) b^{l'}_\alpha, Q_\beta(\theta_k^{l},\theta^l,\cdot) b_\beta^{l} \rangle_\mathcal{H}, \end{align}\] where \(b^l_1 = b(\theta_k^l,\cdot)\) and \(b^l_2=b(\theta^l,\cdot)\). These inner products are then computed either analytically (via closed form formulas) or numerically, depending on the choice of basis functions and the functions \(\phi\) and \(\vec{v}\). For the components of the gradient \[\nabla_{\theta} F_k(\theta) = [ \nabla_{\theta^1} F_k(\theta), \ldots, \nabla_{\theta^L} F_k(\theta) ],\] the above summations shrink, and we have \[\begin{align} \nabla_{\theta^r} F_k(\theta) = \sum_{l'=1}^L & \left( \nabla_{\theta^r} \langle Q_1(\theta^{l'}_k,\theta^{l'},\cdot) b(\theta^{l'}_k,\cdot), Q_2(\theta^r_k,\theta^r,\cdot) b(\theta^r,\cdot) \rangle_\mathcal{H}\right. \\ & + \left. \nabla_{\theta^r} \langle Q_2(\theta^{l'}_k,\theta^{l'},\cdot) b(\theta^{l'},\cdot), Q_2(\theta^{r}_k,\theta^{r},\cdot) b(\theta^r,\cdot) \rangle_\mathcal{H}\right), \quad {\rm for} r=1,\ldots,L. \end{align}\]

3.3.4 Gaussian basis functions↩︎

If a Gaussian basis is used and \(\mathcal{H}=L^2(\mathbb{R}^d,\mathbb{C})\), then the cost function \(F_k(\theta)\) and its gradient \(\nabla F_k(\theta)\) can be computed in two different ways, depending on whether the function \(\phi\) in 9 is polynomial or not.

If \(\phi\) is a polynomial, then the integrals defining \(F_k(\theta)\) are computed explicitly via polynomial Gaussian integrals. More precisely, given any polynomial \(\phi\), the integral \(F_k(\theta)\) and its gradient \(\nabla_\theta F_k(\theta)\) are computed in terms of sums of monomial Gaussian integrals, see 5. An important advantage of Gaussian approximations is that the computation of these integrals in closed form can be generalized to higher dimensions. With analytic formulas, the cost of residual evaluation scales quadratically with the number \(L\) of basis functions. The dominant computational expense is then the nonlinear optimization.

If \(\phi\) is not a polynomial, then we compute each summand \(\langle Q_\alpha b^{l'}_\alpha, Q_\beta b_\beta^{l} \rangle_{L^2}\) via quadratures (we omit the \(\theta\)-dependence for readability), i.e., \[\begin{align} \langle Q_\alpha b^{l'}_\alpha, Q_\beta b_\beta^{l} \rangle_{L^2} &\approx \sum_{x_j\in {\rm Grid}} \overline{Q_\alpha(x_j) b^{l'}_\alpha (x_j)} Q_\beta (x_j) b_\beta^{l} (x_j) w_j, \\ \nabla_{\theta^{l}} \langle Q_\alpha b^{l'}_\alpha, Q_\beta b_\beta^{l} \rangle_{L^2} &\approx \sum_{x_j\in {\rm Grid}} \overline{Q_\alpha(x_j) b^{l'}_\alpha (x_j)} \nabla_{\theta^{l}} \left[ Q_\beta (x_j) b_\beta^{l} (x_j) \right] w_j, \end{align}\] for weights \(w_j\ge0\) and grid points \(x_j\in\mathbb{R}^d\). We note that for the gradient formula a holomorphic parametrization of the Gaussians has to be used. For instance, we may consider Gauss–Hermite quadrature if \(d=1\), see 5.2. In higher dimension \(d\ge 1\), tensor product quadrature comes with a computational complexity that increases exponentially with the dimension \(d\); see also [22] for a sparse Gaussian method.

4 Unified Error Analysis↩︎

For the error analysis, we return to the setting of a general Hilbert space \(\mathcal{H}\). For readability, we omit the sub-indices and simply write \(\langle\cdot,\cdot\rangle\) for the inner product and \(\Vert \cdot \Vert\) for the induced norm. For our analysis, the knowledge on whether \(\psi(t,\cdot)\) is approximated by one or several basis functions or by which type of basis functions is not needed. For the general time marching scheme \[\begin{align} \label{eq:theta95k431} \theta_{k+1} &= \arg \min_{\theta\in\Theta} \Vert r_k(\theta,\cdot) \Vert, \\\nonumber u_{k+1} &= u(\theta_{k+1},\cdot),\quad k=0,1,\ldots, \end{align}\tag{12}\] we consider the two types of residual functions in the form \[\label{eq:residual} r_k(\theta,\cdot) = \left\{ \begin{array}{l} u(\theta,\cdot) - u(\theta_k,\cdot) - h \left[ \zeta f(t_k,\cdot,u(\theta_k,\cdot)) + \hat{\zeta}f(t,\cdot,u(\theta,\cdot)) \right],\\*[1ex] \left[ \zeta \nabla_{\theta} u(\theta_k,\cdot) + \hat{\zeta} \nabla_{\theta} u(\theta,\cdot) \right] (\theta - \theta_k) - h \left[ \zeta f(t_k,\cdot,u(\theta_k,\cdot)) + \hat{\zeta} f(t,\cdot,u(\theta,\cdot)) \right], \end{array}\right.\tag{13}\] with parameters \(\theta,\theta_k\in\Theta\) and \(\zeta\in[0,1]\) and \(\hat{\zeta}=1-\zeta\). The two residuals – one of them gradient-free, the other one containing the gradient of \(u\) – correspond to the discretize then parametrize approach 5 and the opposite one 7 , respectively. The residual \(r_k(\theta,\cdot)\) enters the error estimates in 1 and 2 via the maximal norm value \[\rho = \max\left\{\Vert r_k(\theta_{k+1},\cdot) \Vert: 0\le t_{k+1}\le T\right\}.\] The second main error contribution is related to the time-regularity of the solution \(\psi(t,\cdot)\) and the choice of \(\zeta\). We will use constants \(c_1,c_2>0\) defined by \[\label{def:const} \begin{array}{ll} c_1 = \left( |\tfrac12-\zeta| + \tfrac14\right) \max_{t\in[0,T]}\|\psi''(t,\cdot)\|, & \text{if}\;\zeta\in[0,1]\;\text{and}\;\psi\in C^2([0,T],\mathcal{H}),\\*[2ex] c_2 = \frac{1}{12} \max_{t\in[0,T]}\|\psi'''(t,\cdot)\|, & \text{if}\;\zeta = \frac{1}{2}\;\text{and}\;\psi\in C^3([0,T],\mathcal{H}). \end{array}\tag{14}\] Given \(\rho\) and \(c_1,c_2\), we will be able to bound the error of the discretize-then approximation schemes, see 4.1. For the parametrize-then approach, see 4.2, we will require additional lower bounds on the smallest singular value of the parametrizing gradient.

4.1 Error analysis for the discretization-then approach↩︎

The error analysis here extends the results of [15], that consider the explicit and implicit Euler method for Lipschitz continuous right-hand side. [19] also uses the discretize-then setting and performs the error analysis for regularized Radau IIA and Gaußimplicit Runge–Kutta methods. Here, we analyze the full family of \(\zeta\)-methods, \(\zeta\in[0,1]\), and show the expected gain in accuracy for \(\zeta=\frac{1}{2}\) provided that the original dynamics are regular enough. Also in the one-sided Lipschitz and the dissipative regime, it is the implicit Euler method and the ones with \(\zeta\le\frac{1}{2}\) that provide the expected estimates. A common feature of currently available error estimates for parametric approximation methods is the presence of residual terms.

4.2 Error analysis for the parametrize-then approach↩︎

Error and stability analysis for the parametrize-then approach has been carried out in [15] for time-continuous dynamics. The analysis of fully discretized, regularized Runge–Kutta schemes in [18] comprises the explicit and implicit Euler method. There, the regularization induces a step-size restriction that substitutes the upper bound ?? and the lower bound on the parametrizing gradient’s singular values ?? . These two bounds are needed here to control the distance of consecutive parameter values. They are the key ingredients for the new error constant \(c_s'\) that distinguishes the error bound of 2 from the Lipschitz case in 1.

4.3 Norm and energy error↩︎

For the time-dependent Schrödinger equation \({\rm i}\partial_t \psi(t) = \hat{H}\psi(t)\), \(\psi(0) = \psi_0\) with a self-adjoint Hamiltonian operator \(\hat{H}\) and initial condition \(\psi_0\in\mathcal{H}\), there are two conserved quantities. The norm and the energy satisfy \[\Vert \psi(t,\cdot)\Vert = \Vert\psi_0\Vert,\qquad \langle \psi(t,\cdot),H\psi(t,\cdot)\rangle = \langle \psi_0,H\psi_0\rangle\] for all \(t\in\mathbb{R}\). The time-continuous Dirac–Frenkel approximation inherits this conservation property, see e.g. [4] or [13]. Here, however, with the residual-minimization schemes considered above, these conservation errors are controlled by the norm of the global error \(\varepsilon_k = \psi(t_k,\cdot) - u(\theta_k,\cdot)\).

5 Gaussian integrals and initialization↩︎

For a Gaussian basis and \(\mathcal{H}=L^2(\mathbb{R}^d,\mathbb{C})\), the cost function \(F_k(\theta)=\Vert r_k(\theta,\cdot) \Vert^2_\mathcal{H}\) for parameter \(\theta\in\Theta\) can be expressed in terms of sums of Gaussian integrals. They are at the core of successful computational methods like the variational multi-configuration Gaussian method (vMCG) [2]. Here, for a simple and conceptual exposition, we focus on the univariate case \(d=1\). Let us recall that we use complex parameters \(\theta=(a,\kappa,\gamma)\in\Theta =\mathbb{H}\times \mathbb{C}^2 \subset \mathbb{C}^3\) and denote the single Gaussian function by \[\label{eq:gaussian95def951d} g(\theta,x) := e^{-ax^2+\kappa x +\gamma},\quad x\in\mathbb{R}.\tag{17}\] Notice that the identity \(g(\theta_1,x) g(\theta_2,x) = g(\theta_1+\theta_2,x)\) holds for any \(\theta_1,\theta_2\in\Theta\) with this parametrization choice, which leads to compact formulas for the inner products.

5.1 Monomial integrals↩︎

For \(n\in\mathbb{N}\), we then define the \(n\)-monomial Gaussian integral as \[\label{eq:Gn} G_n(\theta) := \int_{\mathbb{R}} x^n g(\theta,x) {\rm d}x.\tag{18}\] The case \(n=0\) is the usual Gaussian integral \[\label{eq:G0} G_0(\theta) = e^{\frac{\kappa^2}{4a}+\gamma} \sqrt{\frac{\pi}{a}},\tag{19}\] which holds for complex parameters \(\theta\) as long as \(a\in\mathbb{H}\), see [23].

With this formula, one can determine the explicit form of the integral \(G_n(\theta)\) for any \(n\in\mathbb{N}\), using complex scaled¹ (physicists’) Hermite polynomials \(\tilde{H}_n(x):=(-{\rm i})^n H_n({\rm i}x)\). They obey the three-term recurrence \[\begin{align} &\tilde{H}_{n+1}(x)=2x \tilde{H}_n(x) + 2n \tilde{H}_{n-1}(x),\quad n>0,\\ &\tilde{H}_0(x)=1, \;\tilde{H}_1(x)=2x. \end{align}\] This scaling provides a convenient representation of the integrals without additional powers of \({\rm i}\).

5.2 Polynomial integrals and quadrature↩︎

For a Gaussian integral with respect to a function \(\varphi:\mathbb{R}\to\mathbb{C}\) we denote \[I(\theta;\varphi) := \int_{\mathbb{R}} \varphi(x) g(\theta;x) {\rm d}x.\] For a vector field \(\vec{\Phi}:\mathbb{R}\to\mathbb{C}^m\) with components \(\Phi_1,\ldots,\Phi_m\), we set \[I(\theta;\vec{\Phi}) := \begin{pmatrix} I(\theta;\Phi_1) & \ldots & I(\theta;\Phi_m) \end{pmatrix}.\] With this compact notation, Gaussian inner products and their parameter gradients can be stated as follows.

For polynomial functions, we obtain a compact description of the Gaussian weighted inner products in terms of the monomial integrals \(G_n(\theta)\), \(n\ge0\), introduced in 18 . For a multi-variate extension of this type of formula, see [24].

When closed-form expressions are unavailable, the integrals must be approximated numerically. For univariate Gaussian integrals with smooth integrands, Gauss–Hermite quadrature \[\label{eq:GH95quadrature} \mathcal{Q}_n(\varphi) = \sum_{j=1}^n w_j \varphi(x_j)\tag{20}\] seems to be the method of choice. In the numerical experiments of 7.4, we have considered integrand functions \(\varphi\), that allow for a change of variables in the complex plane, and have used the approximation \[I(\theta;\varphi) \approx \mathcal{Q}_n (\theta;\varphi)\] with \[\label{eq:GH95quad95theta} \mathcal{Q}_n (\theta;\varphi) = \frac{e^{\kappa^2/4a+\gamma}}{\sqrt{a}} \sum_{j=1}^n w_j \varphi\left( \frac{x_j}{\sqrt{a}} + \frac{\kappa}{2a} \right).\tag{21}\]

5.3 Initializing a Gaussian mixture↩︎

The initial parameter vector \(\theta_0\in\Theta\) determines the initial approximation \(u(\theta_0,\cdot)\in\mathcal{M}\) of the given datum \(\psi_0\). As highlighted by the error estimates of Theorems 1 and 2, the quantity \[\|\psi_0-u(\theta_0,\cdot)\|_{\mathcal{H}}^2\] constitutes one of the relevant error sources. A sufficiently accurate initialization is therefore essential for observing the asymptotic behavior predicted by the theory.

5.3.0.1 Gaussian notation.

In this section, for a clear geometric interpretation, we will denote the Gaussians in terms of their width coefficient \(a\in\mathbb{H}\), center \(\nu\in\mathbb{C}\) and amplitude \(C\in\mathbb{C}^*\) and to avoid confusion with the previous use of \(\theta\), we denote this parametrization by \(\xi\). Namely, we denote the Gaussian basis by \[b(\xi,x) = C e^{-a(x-\nu)^2},\quad x\in\mathbb{R},\] with \(\xi=(a,\nu,C)\). The theta parameter from the previous subsection in terms of these data is \(\theta=(a,2a\nu,-a\nu^2+\log C )\).

For the initial Gaussian approximations \(u(\xi_0,\cdot) = \sum_{l=1}^L b(\xi_0^l,\cdot)\) with \(\xi_0^l = (a_0^l,\nu_0^l,C_0^l)\) used in the numerical experiments of Sections 6 and 7, we have employed the heuristic construction described below. It assumes that the critical points of \(\psi_0\) are known. The cost function under consideration in the following construction is \[F(\xi) = \|\psi_0-u(\xi,\cdot)\|_{\mathcal{H}}^2,\] with the parameter space \(\Xi=(\mathbb{H}\times\mathbb{C}\times\mathbb{C}^*)^L\).

5.3.1 Basic initialization↩︎

The objective of the following procedure is not to solve the global optimization problem \(\xi_0=\arg\min_{\xi\in\Xi} F(\xi)\), which may be ill-posed or characterized by many local minima, but rather to generate a well-conditioned initial guess with subsequent local optimization. In addition to approximation quality, the construction seeks to avoid strong overlap between basis functions, poor conditioning and thus adversely affects on the parametrize-then formulation.

• Select Gaussian number and width. Set the number of basis functions \(L\geq 1\) and set the width parameters \(a_0^l=a\) for all \(l=1,\ldots,L\) for some constant \(a>0\). Of course, the width need not be uniform and could be chosen in temperate progression as in [25]. The number of basis functions \(L\) could be the number of critical points of \(\psi_0\).

• Place Gaussian centers. Center the Gaussians at the critical points of \(\psi_0(x)\), which we denote as \(\nu_0^l\) for \(l=1,\ldots,L\).

• Determine amplitudes. Choose \(\gamma_0^l\) as the solutions of the set of equations \(\psi_0(\nu_0^l) - u(\xi,\nu_0^l)=0\) with \(l=1,\ldots,L\).

• Local optimization. Improve the approximation \(\xi_0\) by minimizing the cost function \(F(\xi)\), using the values of steps (i)–(iii) as the initial guesses, i.e., \(\xi^{(0)}=(a_0^l,\nu_0^l,C_0^l)_{l=1,\ldots,L}\).

For quantum mechanical problems, we assume that \(\Vert\psi_0\Vert=1\) and the local optimization problem will be constrained to \(\Vert u(\xi,\cdot)\Vert=1\).

5.3.2 Refinement strategy↩︎

If the attained residual remains above a prescribed tolerance, repeat steps (i)–(iv) with a larger number of Gaussian basis functions. The refinement is performed by replacing each Gaussian component by an odd number \(m\ge 1\) of localized sub-components whose centers are distributed according to the roots of an appropriate Hermite polynomial. Since derivatives of Gaussians are Gaussian-envelope polynomials, these roots provide natural locations for the refinement, see also 5.3.3. Specifically, for a Gaussian component \(b(\xi_0^l,\cdot)\) with \(\xi_0^l=(a_0^l,\nu_0^l,C_0^l)\), the refined centers are initialized as \[\nu_0^{l'} = \nu_0^l+\frac{x_{l'}}{\sqrt{a_0^l}}, \qquad l'=1,\ldots,m,\] where \(x_{l'}\) are the roots of the physicists’ Hermite polynomial \(H_m\). The remaining parameters are then obtained from the optimization problem of step (iv), subject to fixed centers and width parameters \(a_0^{l'}\) subject to the linear constraint \((1+\varepsilon)a^l_0 \leq a^{l'} \leq a_{\max}\) for some small \(\varepsilon>0\) and an upper bound \(a_{\max}>(1+\varepsilon)a_0^l\). Without the constraints, the optimizer will yield a single Gaussian as the trivial solution, e.g., the centers of the sub-components will cluster around the same point or one sub-component will converge to the original Gaussian and the rest will take weight zero. For all the experiments presented in this paper, we have used initial constants \(a=1.5\) for step (i) of the basic initialization and initial widths parameters \(a_0^{l'}=1.5 a_0^l\) for the Gaussian subcomponents (refinement stage), which we also use as the upper bound, i.e., \(a_{\max}=1.5 a_0^l\).

5.3.3 Polynomials with a Gaussian envelope↩︎

The initialization described above relies on the availability of characteristic locations of the target function, most notably its critical points. For general initial data these points may have to be determined numerically. For the class of Gaussian-envelope polynomials, however, they can be computed systematically. This class is relevant for the numerical examples considered below and is preserved under differentiation. We consider initial data of the form \[\psi_0(x)=P_n(x)b(\xi^{\rm e},x), \qquad x\in\mathbb{R},\] where \(P_n\) is a polynomial with complex coefficients and \(b(\xi^{\rm e},\cdot)\) a Gaussian with parameters³ \(\xi^{\rm e}=(a^{\rm e},\nu^{\rm e},C^{\rm e})\in \mathbb{H}\times\mathbb{C}\times \mathbb{C}^*\). Let \(P_m\), \(m\ge n\), be defined recursively by \[P_{m+1}(x) = 2a^{\rm e}(x-\nu^{\rm e})P_m(x)-P_m'(x), \qquad m\ge n,\] with initial polynomial \(P_n\) given by the representation of \(\psi_0\). Then repeated differentiation gives \[\psi_0^{(j)}(x) = (-1)^j P_{n+j}(x)b(\xi^{\rm e},x), \qquad j\ge 0 .\] Since the Gaussian factor does not vanish, the zeros of \(\psi_0^{(j)}\) coincide with the zeros of \(P_{n+j}\).

For initial data of this form, the locations required in the initialization procedure of 5.3.1 are therefore obtained from the roots of \(P_{n+1}\). Likewise, the roots needed in the refinement strategy of 5.3.2 arise from the same recursion applied to individual Gaussian components. In the special case \(P_n\equiv 1\), the recursion generates scaled Hermite polynomials. Thus, the Hermite-root placement used in the refinement step can be viewed as the pure-Gaussian instance of the same general mechanism.

5.3.4 Illustration of the initialization procedure↩︎

We illustrate the initialization strategy for two of the initial conditions of the numerical experiments in 6 and 7. The aim is to construct a compact Gaussian mixture that accurately represents the initial condition while providing sufficient expressiveness for the subsequent time evolution. Since the dynamics may rapidly leave the single-Gaussian manifold, the initialization is guided not only by the approximation quality at time \(t=0\) but also by the anticipated complexity of the evolution.

2 illustrates the initialization procedure. Here, we show the approximation with three Gaussian components to emphasize the improvement with step (iv), since using more Gaussians already provides a very good approximation with only steps (i)–(iii). The critical-point construction provides suitable locations for the Gaussian components, and the subsequent constrained optimization improves the approximation without introducing degenerate Gaussian configurations. The resulting mixtures serve as initial values for the numerical experiments in Sections 6 and 7.

6 Qualitative Schrödinger dynamics↩︎

We conduct numerical experiments for the uni-variate Schrödinger equations \[{\rm i}\partial_t\psi(t,x) = \left( -\frac{1}{2}\partial_x^2 + V(t,x)\right) \psi(t,x),\quad \psi(0,x) = \psi_0(x)\] with three choices for the potential function \(V(t,x)\) and various initial conditions \(\psi_0\in \mathcal{H}\) with \(\mathcal{H}= L^2(\mathbb{R},\mathbb{C})\). The purpose of the present section is to illustrate the qualitative behavior for several representative Schrödinger dynamics. Quantitative investigations of residual minimization, time discretization, quadrature errors, norm and energy conservation are presented in 7.

6.0.0.1 Experimental design.

The error estimates of [prop:discr-opt,prop:opt-discr] separate the contributions of time discretization and residual minimization. In practical computations, however, additional sources of error are present. These include the initialization error as well as quadrature errors arising when residual norms and gradients are evaluated numerically. Throughout the experiments, these additional errors are controlled sufficiently so that the asymptotic behavior predicted by the theory becomes observable. Furthermore, the number \(L\) of Gaussian basis functions was chosen such that the approximation manifold is sufficiently expressive while maintaining a reasonably well-conditioned parametrization. Since the parametrize-then-discretize formulation depends explicitly on the conditioning of the parametrization gradient, this choice enables a meaningful comparison of the two residual formulations under conditions where neither method is a priori disadvantaged.

We have considered the following three test cases with different initial conditions.

6.0.0.2 Test 1: Harmonic oscillator (exactly representable benchmark).

The potential is the standard harmonic one, \[V(t,x) = \frac{1}{2} x^2,\quad x\in\mathbb{R},\quad t\in[0,10].\] We recall that the ground state is the centered Gaussian \(\psi_{\rm GS}(x)=\mathcal{N}_0 e^{-\frac{1}{2} x^2}\) with energy \(E_{\rm GS}=\langle\psi_{\rm GS},\hat{H}\psi_{\rm GS}\rangle_\mathcal{H}=\frac{1}{2}\). We consider two different initial conditions, for which harmonic quantum motion can be explicitly resolved within a nonlinear parametrization. The first one lies within the range of a Gaussian parametrization, whereas the second one requires Hagedorn functions, see [3] or the review [26].

Shifted Gaussian.

The shifted Gaussian \(\psi_0(x)=\mathcal{N}_0 e^{-(x-2)^2}\) has its energy \(E_0 = \langle\psi_0,\hat{H}\psi_0\rangle_\mathcal{H}= 2.625\) significantly above the ground state energy and has a slimmer initial width than that of the harmonic ground state. 3 shows the evolution of this state. The numerical solution indicates that the width coefficient \(\Re(a)\) oscillates in the interval \([\frac{1}{4},1]\) and its center oscillates (with a different period) in the interval \([-2,2]\).

Shifted excited state.

The initial condition \(\psi_0(x)=\mathcal{N}_0 H_2(x-X_0) e^{-\frac{1}{2}(x-X_0)^2}\) is similar to the second excited state, except for the initial shift to the right given by \(X_0=1\). We evolve this state considering a three-Gaussian approximation, see 5.3.3 for the determination of the initial value \(\theta_0\). 3 shows that the constant-time sections of the density \(|\psi(t,x)|^2\) remains identical to that of the first excited state \(\psi_{1S}\), except for a shift by \(X\in [-X_0,X_0]\); namely, the probability density is of the form \(|\psi(t,x)|^2=|\psi_{1S}(0,x-X(t))|^2\).

This example serves as a benchmark case in which the Gaussian manifold is invariant under the exact dynamics. Consequently, the residual can be minimized to machine precision and the observed errors are dominated by time discretization.

6.0.0.3 Test 2: Double-well (interference and splitting).

The potential is a double-well \[V(t,x) = \alpha_4 x^4 + \alpha_2 x^2,\quad x\in\mathbb{R},\quad t\in[0,10].\] Following [18], [27], we choose the values for the potential as \(\alpha_2=-1/8\) and \(\alpha_4=(\alpha_2)^2\). The local minima \(y_{\min}=-1/4\) are then located at \(x_{\pm}=\pm 2\), and the local maximum \(y_{\max}=0\) is located at \(x_0=0\). We note that double-well dynamics cannot be exactly resolved within nonlinear parametrization.

First excited state.

We start by studying the evolution of the first excited state in order to visualize the stability of the numerical schemes presented in the paper. To approximate this state, we consider the two-Gaussian approximation \(u_0(x) = \mathcal{N}_0 ( e^{-a_*(x-\nu_*)^2} - e^{-a_0(x+\nu_*)^2} )\), where \(a_*=0.4475\) and \(\nu_*=1.8732\) were determined numerically by minimizing the energy function \(E(\theta)=\langle u(\theta,\cdot),\hat{H} u(\theta,\cdot) \rangle\) and imposing an antisymmetry condition on the Gaussian components of \(u(\theta,\cdot)\). The energy of this state is \(E_{0} = 0.1668\). The numerical result shown on the left-hand side of 4 shows the stationary nature of this state, i.e., the probability density function is time-independent.

Gaussian superposition.

Then we continue with the superposition of Gaussians centered in the left and the right well, namely, \(\psi_0 (x) = \mathcal{N}_0 ( e^{-(x-2)^2} - e^{-(x+2)^2} )\). The energy of this state is \(E_0=0.3181\), which is larger than the first excitation energy. The centers of the Gaussian components are shifted relative to the centers \(\pm\nu_*\) of the first excited state, and their initial width is significantly slimmer.

This example illustrates a regime in which the solution leaves the single-Gaussian manifold and motivates the use of multiple Gaussian components.

6.0.0.4 Test 3: Hyperbolic cosine (strongly non-polynomial dynamics).

The potential is the hyperbolic cosine \[V(t,x) = \cosh(x)-1, \quad x\in\mathbb{R}, \quad t\in [0,10].\] We choose this potential since its second-order Taylor expansion about \(x=0\) is the harmonic oscillator. However, the dynamics will be different in general. We study the evolution of a single Gaussian initial condition \(\psi_0(x)=\mathcal{N}_0 e^{-\frac{1}{2}(x-1)^2}\) considering a multi-Gaussian approximation, i.e., we consider a refinement as described in 5.3.2. For reference, the energy of this state is \(E_0=1.2313\), which is larger than the energy \(E_{GS}=0.5301\) of the single-Gaussian approximation of the ground state \(\psi_{GS}(x)= \mathcal{N}_0 e^{-0.5591 x^2}\). The same initial condition in a harmonic oscillator potential would produce a solution whose density \(|\psi(t,x)|^2\) constant-time sections do not change, except for a shift in the center of the Gaussians \(X\in [-1,1]\), as in the example shown in the left-hand side of 3. This is no longer the case in a hyperbolic cosine potential, as one can observe in 5. Moreover, one single Gaussian is not enough to describe the correct behavior of the wave function, but a refinement with at least three Gaussians is needed. The plot on the right panel shows the evolution for a \(L=5\) approximation.

This example demonstrates that the framework remains effective even when the dynamics are no longer generated by a polynomial Hamiltonian and exact representability is lost.

6.0.0.5 Reproducibility: algorithmic parameters.

All numerical experiments have been done using a unified Matlab code⁴ admitting the following choices: (i) the discretization/parametrization approach, (ii) the weight \(\zeta\in [0,1]\), (iii) the step-size \(h>0\), (iv) the computation of the inner products using quadratures or closed-form formulas (if the potential is a polynomial). The optimization part described in 2.3 has been done with Matlab’s fminunc built-in function with a fixed maximum number of iterations MaxIterations=\(200\) and a constant value \(10^{-8}\) for the tolerances OptimalityTolerance, StepTolerance, FunctionTolerance. Unless explicitly stated, all the experiments shown in this section have been done using the following choices: (i) discretize-then approach, (ii) weight \(\zeta=\frac{1}{2}\), (iii) step-size \(h=10^{-2}\), (iv) closed-form formulas for the polynomial potentials and quadratures for the non-polynomial one.

7 Quantitative numerical validation↩︎

The purpose of this section is to validate the predictions of Theorems 1 and 2. as well as Proposition 3. In particular, we investigate the dependence of the global approximation error on the residual magnitude, the time-step size \(h\), the parameter \(\zeta\), and the accuracy of the quadrature approximation used in the residual evaluation. As quantitative measures we will use the value of the cost function at the optimized parameter \(\theta_k\) for time step \(t_k = kh\), \(k\ge 0\), as well as the norm and the energy deviation, \[F_{k}(\theta_{k+1}) = \Vert r_{k}(\theta_{k+1},\cdot)\Vert^2,\quad \left|\Vert u_k\Vert -1\right|,\quad \left|E_k/E_0-1\right|,\] where \(u_k = u(\theta_k,\cdot)\) and \(E_k = \langle u_k,\hat{H}u_k\rangle_\mathcal{H}\) is the energy associated with the Hamiltonian \(\hat{H}=-\frac{1}{2}\partial_x^2 + V(x)\). We consider maximal and averaged quantities, e.g. \[\max_{k\in K} F_{k}(\theta_{k+1}),\quad {\rm avg}_{k\in K} F_k(\theta_{k+1}) = \frac{1}{N} \sum_{k=1}^N F_{k}(\theta_{k+1})\] with \(N=\lceil T/h \rceil\) the number of time steps \(t_k\) until final time \(T\) and \(K = \{1,\ldots,N\}\). We will see that maximal and averaged quantities are comparable for all our experiments.

7.1 Influence of the manifold dimension↩︎

Increasing the number of Gaussian basis functions enlarges the approximation manifold and therefore reduces the residual contribution appearing in the error estimates of Theorems 1 and 2. 6 illustrates this effect for two of the test cases: the double-well with initial condition \(\psi_0(x)=\mathcal{N}_0 e^{- (x-2)^2}\) and the hyperbolic cosine potential with \(\psi_0(x)=\mathcal{N}_0 e^{-\frac{1}{2} (x-1)^2}\). While a single Gaussian already captures the dominant localization of the wave packet, a small number of additional components improves the approximation of its finer structure.

7.2 Influence of the parameter \(\zeta\)↩︎

The exact representability of the harmonic oscillator suggests that the midpoint choice \(\zeta=\frac{1}{2}\) is optimal whenever the residual vanishes. This situation occurs for the harmonic oscillator, where the Gaussian manifold is invariant under the exact dynamics. Numerical experiments readily confirm this prediction (results not shown). We now investigate whether the midpoint choice remains advantageous when the residual cannot be eliminated completely, which is the generic situation. For this, we consider the hyperbolic cosine and the double-well potential. For both experimental set-ups, the parameter \(\zeta\) is sampled on a grid in \([0,1]\) with increased resolution near \(\zeta=\frac{1}{2}\). The final time is \(T=2\) and the time step-size is fixed at \(h=10^{-2}\).

7.2.0.1 Hyperbolic cosine potential.

The results for the hyperbolic cosine potential use the initial condition \(\psi_0(x)=\mathcal{N}_0 e^{-\frac{1}{2} (x-1)^2}\) and work with a three-Gaussian approximation, \(L=3\). They are shown in 7. In contrast to the harmonic oscillator, the residual cannot be driven to zero; the (averaged) attained minima of the cost function are of order \(10^{-7}\) for both residual formulations. However, we note that the cases \(\zeta>0.5\) did not converge for the discretize-then formulation. There, a finer step-size than \(h=10^{-2}\) is needed for convergence. Despite these comparable residual levels (for the cases that converge), the resulting norm and energy deviations exhibit a pronounced dependence on \(\zeta\). In particular, the parametrize-then approach produces smaller norm and energy deviations for almost all values of \(\zeta\) (except around \(\zeta=\frac{1}{2}\), where discretize-then is visibly better). For both the discretize-then and the parametrize-then formulations, the smallest deviations are attained for values of \(\zeta\) close to \(\frac{1}{2}\). In the discretize-then formulation, the minimum is achieved at \(\zeta=\frac{1}{2}\), whereas for the parametrize-then formulation it occurs at a nearby value.

7.2.0.2 Double-well potential.

In 8, we repeat the same experiment for the double-well potential and the initial condition \(\psi_0(x)=\mathcal{N}_0 e^{-(x-2)^2}\), which is centered at the bottom of the right well. Since a single Gaussian will not properly resolve the dynamics, we choose \(L=3\) for the approximation. It can be observed that the attained values of the cost function are around \(10^{-6}\) and \(10^{-7}\). As in the preceding example, the discretize-then approach will not converge in the vicinity of \(\zeta=1\) for the chosen step-size \(h=10^{-2}\). A smaller step-size is needed in such cases. We also observe that the optimal choice for \(\zeta\) is again \(\zeta\approx \frac{1}{2}\). However, unlike the hyperbolic cosine example, the parametrize-then formulation does not perform better than discretize-then in all situations away from the midpoint. In this example, the discretize-then approach produces slightly smaller norm deviations for most values of \(\zeta\) (except around \(\zeta=\frac{1}{2}\)), but parametrize-then produces smaller energy deviations.

7.2.0.3 Conclusion.

In agreement with our theoretical findings, choosing \(\zeta\approx \frac{1}{2}\) is the preferred option for both schemes (discretize-then/parametrize-then), even in situations in which the cost function cannot be minimized to zero. Moreover, the examples above show situations in which either of the schemes has a better performance. In particular, it seems that discretize-then is the optimal choice when \(\zeta=\frac{1}{2}\), whereas parametrize-then is the best choice when \(\zeta=1\).

7.3 Convergence with respect to the step-size↩︎

We next investigate the dependence of the approximation error on the time-step size \(h\). According to Theorems 1 and 2, the global error consists of two contributions: a time-discretization error determined by the chosen \(\zeta\)-scheme and a residual contribution associated with the restriction to the approximation manifold. The purpose of the present experiment is to illustrate how these two contributions interact.

For the harmonic oscillator, the Gaussian manifold is invariant under the exact dynamics. Consequently, the residual can be reduced to machine precision, provided that the step-size is sufficiently small, and the observed error is governed essentially by the time discretization. In this residual-free regime, the expected orders are visible: Euler-type choices of \(\zeta\) lead to first-order convergence, whereas the midpoint choice \(\zeta=\frac{1}{2}\) yields second-order convergence (no numerical results presented here).

The situation changes when the exact dynamics are no longer contained in the Gaussian manifold. To illustrate this effect, we consider the hyperbolic cosine potential, for which the residual contribution in the error estimates does not vanish. Figure 9 shows the scaled residual \(\rho/h = \max_{k\in K} \sqrt{ F_k(\theta_{k+1})}/h\) and the maximal norm deviation \(\Delta u_{\max}:= \max_{k\in K} \bigl|\Vert u(\theta_k,\cdot)\Vert-1\bigr|\) as a function of the step-size \(h\) for the discretize-then and parametrize-then formulations. The computations use step-sizes \(h\in[10^{-2},10^{-1}]\) with equidistant values of \(\log h\), final time \(T=2\), and Gaussian mixtures with \(L=1,3,5,7\) components. Table 2 reports the corresponding regression parameters for the models \[\Delta u_{\max} \approx C_u h^{p_u}.\] The results show the expected first-order behavior for the implicit Euler case \(\zeta=0\). The fitted exponents are close to one in most cases, indicating that the contribution of the time discretization decreases linearly with \(h\). For \(\zeta=\frac{1}{2}\), the interpretation is more delicate. The midpoint rule reduces the discretization error, but the residual contribution does not vanish. Once this residual contribution becomes comparable to, or larger than, the time-discretization error, the second-order regime is partially masked. This explains why the observed slopes do not always display a clean quadratic rate, even though the midpoint rule still produces substantially smaller norm deviations.

This behavior is particularly visible in the discretize-then formulation, where for \(\zeta=\frac{1}{2}\) the norm deviation reaches its minimal level already for relatively large step-sizes. Further decreasing \(h\) then produces little improvement, because the dominant error is no longer the time-discretization error but the residual contribution associated with the Gaussian manifold. In contrast, the parametrize-then formulation shows a clearer decrease in some cases, most notably for \(L=7\), where the fitted exponent is close to two. This indicates that, when the residual is sufficiently well controlled, the second-order behavior of the midpoint rule can become visible.

7.3.0.1 Conclusion.

Overall, Figure 9 and Table 2 are consistent with the error decomposition predicted by the theory. A reduction of the step-size improves the approximation only until the residual contribution becomes dominant. The midpoint choice \(\zeta=\frac{1}{2}\) remains advantageous, not necessarily because it always improves the observed asymptotic order, but because it consistently reduces the error constants and yields smaller norm deviations over the tested range of step-sizes. Similar observations apply for the double-well dynamics (numerical experiments not presented here).

7.4 Influence of quadrature accuracy↩︎

We compare the obtained results when the residual function is computed with closed-form formulas (only possible if the potential has a special form, for example, polynomial) and those obtained with Gauss–Hermite quadrature as described in 5.2. We choose the rule with \(N=10\) nodes, which is exact for integrands that are a polynomial times a Gaussian for polynomial degree \(\le 2N-1\). However, this assumes that the nodes and weights are exact, but in practice, they can only be obtained numerically. The residual evaluation involves repeated optimization steps and repeated quadrature evaluations. Small inaccuracies are therefore propagated through the nonlinear optimization process and accumulate over time.

In the harmonic oscillator example ( [fig:ex-vs-quad-a,fig:ex-vs-quad-b]), the quadrature and exact evaluation results are almost indistinguishable throughout the simulation interval. In the double-well example ([fig:ex-vs-quad-e,fig:ex-vs-quad-f]), small discrepancies appear rather early and slightly accumulate over time. All double-well results reflect the greater sensitivity of the generic case, when the approximation manifold is not invariant under the dynamics.

These experiments indicate that, in one spatial dimension, a \(10\)-point Gauss–Hermite quadrature is sufficiently accurate to reproduce both the residual evolution and the conservation properties obtained with exact Gaussian formulas. Nevertheless, exact formulas remain preferable whenever available, since they eliminate quadrature error entirely and become increasingly advantageous in higher-dimensional settings.

8 Conclusions↩︎

We have presented a unified analytical framework for time-discrete parametric approximations of evolution equations in Hilbert spaces based on residual minimization for nonlinear approximation manifolds. The formulation encompasses two natural paradigms: discretization followed by parametrization of the evolution equation (Rothe approach), and discretization of the variational dynamics induced by the Dirac–Frenkel principle. Both strategies are expressed in terms of residual norms whose minimization defines the time-stepping scheme in parameter space.

The main theoretical contribution is a unified error analysis for the full family of \(\zeta\)-methods, extending previous results in [15] and complementing those for regularized parametric approximation [18], [19]. The resulting bounds separate the effects of time discretization from those arising by residual minimization. This structure is independent of the particular form of the underlying evolution equation and relies only on general Lipschitz-type or dissipativity properties. For the parametrize-then formulation, additional conditioning requirements on the parametrization appear naturally and explain stability constraints encountered in nonlinear manifold approximations. For the discretize-then approach with implicit Euler, the condition can be relaxed to a one-sided Lipschitz condition. Discretize-then with \(\zeta<\frac{1}{2}\) and a one-sided Lipschitz condition with \(\ell=0\) allows proving analogous bounds.

While Gaussian approximation manifolds provide a concrete and computationally attractive realization of the abstract framework, the analysis itself is not tied to any particular functional form. The residual-based viewpoint applies equally to a broad range of nonlinear parametrizations, including reduced-basis spaces, adaptive dictionaries, and data-driven manifolds. The numerical experiments for time-dependent Schrödinger equations serve primarily as illustrative examples of the general theory. They confirm the predicted convergence orders and highlight the role of residual accuracy. The Fokker–Planck example presented in Appendix  12 further illustrates that the residual-based approximation is not tied to Hamiltonian dynamics and can be applied to dissipative evolution equations as well.

The flexibility of the framework also suggests several directions for future research. In particular, adaptive approximation manifolds appear especially promising from a computational perspective. Rather than working with a fixed approximation space, one may allow the manifold to evolve together with the solution, enriching or simplifying the representation whenever required by the local approximation quality. In the Gaussian setting, this corresponds to dynamically adjusting the number of basis functions during the computation. More generally, the residual-based formulation developed here may provide a natural foundation for the analysis of adaptive parametrizations in a variety of reduced-order and data-driven approximation methods.

Acknowledgements↩︎

Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – TRR 352 – Project-ID 470903074.

Data Availability↩︎

All code used for the qualitative and quantitative numerical experiments is openly available at https://gitlab.lrz.de/00000000014C092A/discretize-then-and-parametrize-then.

9 Local error estimate for \(\zeta\)-methods↩︎

10 Nonlinear least squares estimate↩︎

11 Gaussian integrals↩︎

11.1 Monomial formula↩︎

For proving Formula 4, we compute \[\begin{align} G_n(\theta) &= \int_\mathbb{R}x^n e^{-a x^2 +\kappa x + \gamma } {\rm d}x \;=\; \frac{\partial^n}{\partial \kappa^n} \int_\mathbb{R}e^{-a x^2 +\kappa x + \gamma } {\rm d}x \\ &\stackrel{\eqref{eq:G0}}{=} e^{\gamma } \sqrt{\frac{\pi}{a}} \frac{\partial^n}{\partial \kappa^n} e^{\frac{\kappa^2}{4a} } \;= \; (2\sqrt{a})^{-n} \tilde{H}_n\left( \frac{\kappa}{2\sqrt{a}} \right) \sqrt{\frac{\pi}{a}} e^{\frac{\kappa^2}{4a} +\gamma } , \end{align}\] where we have used the Rodrigues formula for the Hermite polynomials. For the gradient we simply observe \[\begin{align} \nabla_{\theta} G_{n}(\theta) = \left( \begin{array}{c} \partial_a G_n(\theta) \\ \partial_\kappa G_n(\theta) \\ \partial_\gamma G_n(\theta) \end{array} \right)^\top = \left( \begin{array}{c} - \int_{\mathbb{R}} x^{n+2} g(\theta,x) {\rm d}x\\ \int_{\mathbb{R}} x^{n+1} g(\theta,x) {\rm d}x \\ \int_{\mathbb{R}} x^n g(\theta,x) {\rm d}x \end{array} \right)^\top. \end{align}\]

11.2 Gauss–Hermite quadrature↩︎

The nodes \(x_1,\ldots,x_n\) of the Gauss–Hermite quadrature rule 20 are the roots of the \(n\)th Hermite polynomial \(H_n(x)\) and the weights are \[w_j = \frac{2^{n+1} n! \sqrt{\pi}}{[H_{n+1}(x_j)]^2},\quad j=1,\ldots,n.\] The classic error representation for a \(2n\) times continuously differentiable integrand function \(\varphi\) is \[\int_{\mathbb{R}} \varphi(x) e^{-x^2}{\rm d}x = \mathcal{Q}_n(\varphi) + \frac{n!\sqrt{\pi}}{2^n (2n)!} \varphi^{(2n)}(\xi),\] for some \(\xi\in\mathbb{R}\), see e.g. [29].

12 Numerical experiments for Fokker–Planck equations↩︎

The purpose of this appendix is merely to illustrate that the residual formulations developed in the main text are not restricted to Schrödinger equations. We thus consider the Fokker–Planck equation with time-dependent coefficients \[\partial_t \psi (t,x) = D(t) \Delta \psi(t,x) - \nabla (w(x,t) \psi(x,t)),\] where \(D(t)\) is the diffusion coefficient and \(w(x,t)=A(t)x+b(t)\) is the drift coefficient. It is clear that the operator \(f(t,x,\cdot)=D(t) \Delta_x - \nabla_x (w(x,t) \cdot)\) acts on the Gaussian basis linearly, i.e., \(f(t,x,b(\theta,x)) = \phi(\theta,t,x)b(\theta,x)\). For the univariate case \(d=1\), the function \(\phi(\theta,t,x)\) is given by \[\phi(\theta,t,x) = D(t) [ -2a + (-2a x + \kappa)^2 ] - A(t) - w(x,t) (-2a x + \kappa).\]

11 shows the evolution of a single Gaussian wave packet \(\psi_0(x) = \mathcal{N}_0 e^{-(x+1)^2}\) in four different regimes: a) \(A=0\), \(b=0\), \(D=0.02\); b) \(A=0.02\), \(b=0\), \(D=0\); c) \(A=0\), \(b=4/10+\sin(2\pi t)\), \(D=0\); d) \(A=0.01\), \(b=4/10+\sin(2\pi t)\), \(D=0.01\). We recall that b) is the diffusion equation and c) is the transport equation. As in the examples in 6, all the numerical solutions were computed using the discretize-then formulation with the parameters \(\zeta=\frac{1}{2}\) and \(h=10^{-2}\).

In this case, the numerical approximation can be compared with the analytical solution. 12 shows the difference \(\Vert \psi(t_k,\cdot)-u(t_k,\cdot)\Vert^2_{L^2}\) in function of time, corresponding to the example in 11 (d), which will decrease with the step-size until reaching the limit due to round-off errors. For the analytical solution, recall that given a single Gaussian wave packet as the initial condition, the solution will remain a single Gaussian. We write this solution in the form \(\psi(t,x)=e^{-\frac{1}{S}(x-\nu)^2+\tilde{\gamma}}\) with coefficients \(\tilde{\theta} =(S,\nu,\tilde{\gamma})\). Then, \(\psi\) is described by the solutions of the ODE system \[\begin{align} S'(t) = 4D+2AS, \quad \nu'(t) = A\nu + b, \quad \tilde{\gamma}'(t) = -A-2D/S. \end{align}\] Moreover, assume that the coefficient \(A\geq 0\) is constant. Then, the solution of the above ODE system with initial conditions \(\tilde{\theta}_0=(S_0,\nu_0,\tilde{\gamma}_0)\) is \[\begin{align} S(t) &= \left\{ \begin{array}{ll} S_0 + 4 \int_0^t D(t') dt', & \text{if } A=0, \\ e^{2At} \left[ S_0 + 4 \int_0^t D(t') e^{-2At'} dt' \right], & \text{if } A>0. \end{array} \right. \\ \nu(t) &= \left\{ \begin{array}{ll} \nu_0 + \int_0^t b(t') dt', & \text{if } A=0, \\ e^{At} \left[ \nu_0 + \int_0^t b(t') e^{-At'}dt' \right] , & \text{if } A>0, \end{array} \right. \\ \tilde{\gamma}(t) &= \tilde{\gamma}_0 -\frac{1}{2} \ln S(t) , \quad \text{if } A\geq 0, \end{align}\] For instance, considering constant coefficients \(A,D>0\) and the time-dependent function \(b(t)=\alpha+\beta\sin(\omega t)\), the solutions would be \[S(t) = e^{2At} \left[ S_0 + \frac{2D}{A}\left( 1-e^{-2At}\right) \right],\] \[\nu(t) = e^{At} \left[ \nu_0 + \frac{\alpha}{A} (1-e^{-At}) + \frac{\beta\omega}{A^2+\omega^2} \right] - \frac{\beta}{A^2+\omega^2}\left[ \omega\cos(\omega t) + A\sin(\omega t) \right] .\]

References↩︎

---

1. Note that the coefficients of this polynomial are not complex. In fact, its coefficients are the absolute values of the coefficients of the usual Hermite polynomial.↩︎

2. In this paper, we have chosen the branch with \(|{\rm Arg}(\sqrt{a})|<\pi/4\), since we want these formulas to reduce to the usual Gaussian integrals when \(a\) is real, i.e., we want a positive \(\sqrt{a}\) as \(\Im(a)\to 0\).↩︎

3. Here, we use superscripts for the Gaussian envelope to avoid confusion with the Gaussian components of the approximation \(u(\xi,\cdot)\).↩︎

4. https://gitlab.lrz.de/00000000014C092A/discretize-then-and-parametrize-then↩︎