The proposed project will develop and disseminate software that employs a componentwise approach to time-stepping to circumvent the difficulties caused by the phenomenon of stiffness. The new methods will provide a more practical alternative to researchers in many branches of science and engineering who rely on the numerical solution of time-dependent partial differential equations (PDEs). The rapid advancement of computing power in recent years has allowed higher resolution models, which has introduced greater stiffness into systems of ordinary differential equations (ODEs) that arise from discretization of time-dependent PDEs. This presents difficulties for explicit and implicit time-stepping methods, due to the coupling of components of the solution corresponding to low and high frequencies, that change at widely varying speeds.
This goal is to be accomplished by continuing the evolution of Krylov subspace spectral (KSS) methods, which were developed by the PI [54] and have been advanced by the PI over the last several years [9,49,50,51,56,59] These methods feature explicit time-stepping with high-order accuracy, and stability that is characteristic of implicit methods. This ``best-of-both-worlds'' combination is achieved through a componentwise approach, in which each Fourier coefficient of the solution is computed using an approximation of the solution operator that is, in some sense, optimal for that coefficient. Initial work on KSS methods [9,37,49,50,51,52,53,54,55,56,57,58,59,60,69] that includes comparisons to many other time-stepping methods, has yielded promising results, in terms of accuracy, stability, and scalability.
These results encourage their continued development, including enhancement of their accuracy and efficiency, as well as expansion of their applicability. To date, KSS methods have been applied almost exclusively to linear variable-coefficient PDEs on $n$-dimensional boxes, for $n=1,2,3$, with either periodic or homogeneous boundary conditions. The proposed project will feature generalization to KSS methods to nonlinear problems, and to problems on non-rectangular domains.
KSS methods will be extended to nonlinear PDEs by combination with exponential propagation iterative (EPI) methods [79,80] that have shown much promise as effective stiff solvers, compared to existing explicit and implicit time-stepping methods [61]. However, EPI methods currently employ Krylov projection methods [39,40,81] to approximate products of matrix functions and vectors that can require high-dimensional Krylov subspaces, especially at higher spatial resolution. KSS methods, on the other hand, work with Krylov subspaces whose dimension is independent of the spatial resolution. The proposed project will include optimization of KSS methods so that EPI methods, when using KSS methods to approximate these products, will be much more efficient and scalable than in their present form.
KSS methods will also be extended to PDEs on more general domains by combination with Fourier continuation (FC) methods [1,7,63,70] which are applicable to such domains. Specifically, the goal is to combine the spatial discretization of these methods with the approach to high-order time-stepping of KSS methods. As KSS methods are particularly effective for evaluating matrix function-vector products for matrices representing differential operators under periodic boundary conditions at high resolution, FC methods use periodic extensions to achieve high accuracy on general domains, and EPI methods use matrix function-vector products to achieve high-order accuracy for nonlinear PDEs, the combination of the three methods is a logical foundation for a more effective approach to overcoming stiffness.
KSS methods represent the first application to time-dependent PDEs of techniques from the emerging research area of ``matrices, moments and quadrature'': the approximation of bilinear forms involving matrix functions by treating them as Riemann-Stieltjes integrals, and then applying Gaussian quadrature rules that are generated by the Lanczos algorithm. However, as beneficial as KSS methods can be in their own right, the broader impacts of the proposed project on the scientific computing community as a whole will stem from an understanding of the body of work from which they arose.
Because of the widespread importance of simulation techniques as evidenced by the variety of applications, the proposed research effort will be accompanied by educational activities designed to
There are also several applications of these techniques outside of numerical linear algebra. Approximation of bilinear forms is useful for approximation of the scattering amplitude [30], and estimating the trace of the inverse has applications in the study of fractals [74,82] lattice quantum chromodynamics [13,75] and crystals [67,68]. This outreach will encourage researchers to extend the application of these techniques to a wider variety of problems, in which computation of bilinear forms arises naturally.
These activities will be complemented by outreach to high school students and teachers in Mississippi. The PI will travel to area mathematics classrooms and work through exploratory projects with students. They will have hands-on experiences with a variety of methods that produce approximate solutions to problems with which they are familiar, such as nonlinear equations and integrals, and be introduced to application areas. These activities and classroom time with a university researcher will give them an idea of how the mathematics that they are learning is applied in a variety of STEM careers.
Furthermore, for both explicit and implicit methods, these issues are exacerbated as the spatial resolution increases. For explicit methods, the time step must increase linearly (for hyperbolic problems) or quadratically (for parabolic problems) as the mesh spacing decreases, due to the CFL limit. For implicit methods, increased spatial resolution results in an increase in the condition number of the system that must be solved during each time step. Going forward, spatial resolution of models will only continue to increase with improvements in the efficiency of memory and processors, thus creating a situation that is worsening over time. These advances in hardware will foster a demand for, and expectation of, greater use of simulation in real-time (e.g., medicine) and/or high-data (e.g., reservoir simulation) applications that existing time-stepping methods will have difficulty meeting. Therefore, a fresh perspective is needed for the development of time-stepping methods that can more efficiently deal with stiffness.
In many cases, the solution of such a system of ODEs can be expressed in terms of the approximate evaluation of a matrix function such as the exponential. Techniques for computing $f(A){\bf v}$, such as those described in [39,40], have been used for several years to solve systems of ODEs, but can encounter the same difficulties with stiffness. Techniques for computing the bilinear form ${\bf u}^T f(A) {\bf v}$ open the door to rethinking how time-stepping is carried out. Because they allow individual attention to be paid to each component of $f(A){\bf v}$ in some basis, even if those components are coupled, they can be used to construct new time-stepping methods, based on KSS methods, that are more effective for large-scale problems.
A survey of the literature on time-stepping methods shows that recent advances mainly come from two strategies for mitigating the difficulties caused by stiffness: adaptive time-stepping (see, for example, [26,33,42,43,44,45,48,71,84]), in which the global time step used to compute the entire solution is varied according to an error estimate and/or a measure of the smoothness of the solution, and local time-stepping (see, for example, [6,10,16,17,19,22,23,34,35,47,62,78]), in which the time step is varied over the spatial domain so as to satisfy local stability criteria. The guiding principle behind both strategies is to limit the degradation of performance due to stiffness by reducing the time step either locally in time (for adaptive time-stepping) or in space (for local time-stepping) where higher spatial resolution is needed. Also worth noting are multiple time-stepping (MTS) methods [15,21,83], that mitigate the difficulties caused by stiffness by using appropriate time steps for stiff and non-stiff parts of the equation.
By contrast, the approach in the proposed project is to confront stiffness directly by modifying the approximate solution operator in a componentwise manner, so that an increase in spatial resolution does not cause undue additional computational expense, even when using a fixed time step. However, this approach should not be seen as an alternative to the abovementioned strategies. KSS methods can benefit from adaptive time-stepping just like any other time-stepping method through the addition of error estimation; this enhancement will be part of the proposed project. Also, the approach taken by KSS methods of using different approximations for each Fourier component is amenable to local time-stepping, except that the locality is in Fourier space rather than physical space. Finally, MTS methods can benefit from applying KSS methods to the stiff parts of the problem.
For example, consider the problem of computing ${\bf w} = e^{-At}{\bf v}$ for a given symmetric matrix $A$ and vector ${\bf v}$. An approach described in [40] is to apply the Lanczos algorithm to $A$ with initial vector ${\bf v}$ to obtain, at the end of the $j$th iteration, an orthogonal matrix $X_j$ and a tridiagonal matrix $T_j$ such that $X_j^T AX_j = T_j$. Then, we can compute the approximation \begin{equation} {\bf w}_j = \|{\bf v}\|_2 X_j e^{-T_j t} {\bf e}_1, \quad {\bf e}_1 = \left[ \begin{array}{cccc} 1 & 0 & \cdots & 0 \end{array} \right]^T. \label{eqlancapprox} \end{equation} As each column ${\bf x}_k$, $k=1,\ldots,j$, of $X_j$ is of the form ${\bf x}_k = p_{k-1}(A){\bf v}$, where $p_n(A)$ is a polynomial of degree $n$ in $A$, it follows that ${\bf w}_j$ is the product of a polynomial in $A$ of degree $j-1$ and ${\bf v}$.
If the eigenvalues of $A$ are not clustered, which is the case if $A$ arises from a stiff PDE, a good approximation cannot be obtained using a few Lanczos iterations, just as a low-degree polynomial cannot accurately approximate an exponential function on a large interval. This can be alleviated using an outer iteration \begin{equation} {\bf w}_j^{m+1} \approx e^{-A\Delta t}{\bf w}_j^m, \quad m = 0, 1, \ldots, \quad {\bf w}_j^0 = {\bf v}, \label{outerit} \end{equation} for some $\Delta t \ll t$. However, this approach is not practical if $\Delta t$ must be chosen very small. Modifications described in [18,65,81] produce rational approximations that reduce the number of Lanczos iterations, but these methods require solving systems of linear equations in which the matrix is of the form $I-hA$, where $h$ is a parameter. Unless $h$ is chosen very small, these systems may be ill-conditioned. In summary, time-stepping methods that approximate the exponential using a polynomial or rational function tend to have difficulty with stiff systems, such as those that arise from PDEs. The solution of a stiff system includes coupled components that evolve at widely varying rates. As a result, explicit methods for such systems require small time steps, while implicit methods typically require the solution of ill-conditioned systems of linear equations that pose difficulties for direct and iterative solvers alike. The difficulties that many time-stepping methods have with stiffness is associated with the interaction between components of the solution of varying frequency, that vary at different speeds. However, it is not this interaction alone that necessarily precipitates the difficulties typically observed in time-stepping methods applied to stiff problems. A contributing factor is the attempt to treat such disparate components of the solution with the same function, whether polynomial or rational, that approximates $e^{-\lambda\Delta t}$. Unfortunately, such approximation is impractical when $\lambda$ varies within a large interval, which is precisely how the eigenvalues of $A$ behave when the problem (\ref{utAu}) is stiff. Componentwise time-stepping, on the other hand, can avoid these drawbacks, because each Fourier coefficient of the solution is computed using a different approximation of the solution operator. Rather than match terms of the Taylor expansion of the exponential, as in methods of Runge-Kutta type, the exponential function is approximated using an interpolating polynomial with frequency-dependent interpolation points. This increased flexibility can be exploited to obtain the advantage of implicit methods, their stability, without having to solve a large system of equations during each time step. However, considering such an approach raises important questions:
To date, this approach to rapidly estimating quadrature nodes has only been applied to specific classes of differential operators with certain boundary conditions [9,69]. In the proposed project, the goal is to develop an algorithm that can automatically carry out analyses of general linear differential operators, specified in terms of their coefficients and differentiation operators. The algorithm will apply to differential 1-D, 2-D or 3-D operators, with periodic boundary conditions as well as homogeneous Dirichlet, Neumann or Robin boundary conditions. Systems of coupled PDEs can also be accommodated, using an approach described in [53,60]. All of these avenues of generalization have been explored to an extent in [9]; the task for this project is to extrapolate to other differential operators based on patterns observed in previous analyses, at least through three iterations of Lanczos or Arnoldi iteration to produce quadrature nodes for time-stepping methods of up to 5th order.
The success of this approach relies on a determination of what is the low-frequency part of the solution, since different approaches are used to compute low- and high-frequency components [9]. In the proposed project, an algorithm for adaptively determining an appropriate cut-off, based on smoothness of the solution, will be developed. In addition, low-frequency analysis of block Lanczos iteration will be performed, similar to the high-frequency analysis described earlier in this section, so that low-frequency components can be computed using KSS methods rather than Krylov projection, thus reducing Krylov subspace dimension even further.
Women | African-Americans | |||
Category / Year | 2000 | 2012 | 2000 | 2012 |
Mathematics and Statistics | 47.8 | 43.1 | 7.7 | 5.4 |
Computer Science | 28.0 | 18.2 | 9.5 | 10.6 |
All fields | 57.3 | 57.4 | 8.3 | 9.9 |
On the other hand, at USM, approximately 61% of the students are women, while 27% are African-American. As a significant percentage of mathematics majors come from these groups (for example, 80% of seniors graduating in 2009 with bachelors' degrees in mathematics were women), the PI will strongly encourage such students to participate in the project. With the support of requested travel funds, graduate and undergraduate students will be encouraged to participate in conferences such as the AMS Annual Meeting and Southeastern Sectional Meetings, the SIAM Annual Meeting, the International Conference on Spectral and High-Order Methods (ICOSAHOM), and SIAM section meetings. The PI will design special topics courses based on the proposed research, in order to acquaint students with the fundamentals of spectral methods, orthogonal polynomials, and Krylov subspace methods. New courses and modifications of existing courses, as well as their broader impacts on the curriculum, will be assessed through course evaluations completed by students, and exit and alumni surveys of mathematics majors that are currently administered by the department. A learning seminar series based on the project will also be organized, with participating students, the PI, and guest lecturers among the speakers, with travel for visitors supported by a combination of project and department funds.