In many computational science and engineering applications, for example,
fluid flow and interface propagation, it is extremely challenging to design
numerical methods that achieve high order of accuracy. The main theme of my
research is the development of numerical schemes that are high-order both in
space and time, for partial differential equations arise in these applications.

A fundamental drawback of high-order linear multistep methods is their restrictive stability properties compared to RK schemes.
In this work, we develop a new stability theory that characterizes unconditional stability of IMEX multistep methods, and design a new
class of IMEX multistep methods that possesses large regions of unconditional stability.
Specifically, we consider IMEX multistep methods applied to a system of ordinary differential equations $u_t = Au + Bu$, where $Au$ is treated implicitly and $Bu$ is treated explicitly. We are particularly interested in
the problem for which both $A$ and $B$ are stiff. Such splittings can arise, for example, in a pressure Poisson reformulations
of the incompressible Navier-Stokes equations. Unconditional stability allows for time step selection solely
based on accuracy consideration.

The new IMEX multistep methods have tunable parameters that control the size of the unconditional stability region, and can achieve unconditional stability for problems that are unstable with current popular schemes such as semi-implicit backward differentiation formulas (SBDF). Our new stability theory generalizes the existing stability theories, that make use of stability regions and matrix spectra of the problem, to IMEX multistep methods, and provides rigorous explanations for stability limitations in SBDF.

The figures above illustrate the unconditional stability regions for our new IMEX multistep methods from order 2 to 5, with different choice of the parameter $\delta$. Note that $\delta=1$ recovers SBDF methods.

The new IMEX multistep methods have tunable parameters that control the size of the unconditional stability region, and can achieve unconditional stability for problems that are unstable with current popular schemes such as semi-implicit backward differentiation formulas (SBDF). Our new stability theory generalizes the existing stability theories, that make use of stability regions and matrix spectra of the problem, to IMEX multistep methods, and provides rigorous explanations for stability limitations in SBDF.

The figures above illustrate the unconditional stability regions for our new IMEX multistep methods from order 2 to 5, with different choice of the parameter $\delta$. Note that $\delta=1$ recovers SBDF methods.

**Related Publications:**- R. R. Rosales, B. Seibold, D. Shirokoff and D. Zhou.
*Unconditional Stability for multistep IMEX schemes - Theory.*SIAM J. Numer. Anal., Vol. 55, No. 5, 2017, pp. 2336--2360. - R. R. Rosales, B. Seibold, D. Shirokoff and D. Zhou.
*Unconditional Stability for multistep IMEX schemes - Practice.*In preparation.

Despite the advantages of locality in time and good stability properties, RK methods have a key drawback:
in general they are affected by order reduction, meaning that the observed convergence order, when applied to stiff ODEs or
initial boundary value problems (IBVP), is less than the formal order of the scheme.

In this work, we analyze the global errors incurred when applying an implicit RK scheme to a general IBVP $u_t = \mathcal{L}u + f(x, t)$ with time dependent boundary conditions. We develop a theoretical framework that yields a fundamental geometric understanding of the mechanism for order reduction. Geometrically, the order reduction stems from the boundary layers, produced by the fact that the scheme is too accurate near the boundary (the approximation error satisfies a singularly perturbed boundary value problem). The amplitude of the boundary layers is limited by the stage order of the scheme, which in general less than the scheme’s order, thus order reduction.

Two remedies are developed to recover full convergence order:

(1) We derive a new condition, the

(2) The modified boundary conditions (MBC) achieve order reduction avoidance by suppressing the singular behavior in the numerical solutions computed from existing RK schemes. These boundary conditions are extracted via a formal power series expansion of the solution in the absence of the the singular behavior triggered by the conventional approach of enforcing boundary conditions.

The figures above show the error convergence for 3rd order schemes for a 1D heat equation example. Left: 3rd order WSO 1 DIRK scheme with conventional b.c. exhibits order reduction. Middle: The same WSO 1 DIRK scheme with MBC recovers full 3rd order convergence in $u$ and $u_x$. Right: 3rd order WSO 2 DIRK scheme recovers 3rd order convergence in $u$.

In this work, we analyze the global errors incurred when applying an implicit RK scheme to a general IBVP $u_t = \mathcal{L}u + f(x, t)$ with time dependent boundary conditions. We develop a theoretical framework that yields a fundamental geometric understanding of the mechanism for order reduction. Geometrically, the order reduction stems from the boundary layers, produced by the fact that the scheme is too accurate near the boundary (the approximation error satisfies a singularly perturbed boundary value problem). The amplitude of the boundary layers is limited by the stage order of the scheme, which in general less than the scheme’s order, thus order reduction.

Two remedies are developed to recover full convergence order:

(1) We derive a new condition, the

*weak stage order*(WSO), for constructing new RK schemes devoid of order reduction. The novelty of this condition is that it*is*compatible with a diagonally implicit Runge-Kutta (DIRK) structure, in contrast to the high stage order which is impossible to achieve with DIRK schemes.(2) The modified boundary conditions (MBC) achieve order reduction avoidance by suppressing the singular behavior in the numerical solutions computed from existing RK schemes. These boundary conditions are extracted via a formal power series expansion of the solution in the absence of the the singular behavior triggered by the conventional approach of enforcing boundary conditions.

The figures above show the error convergence for 3rd order schemes for a 1D heat equation example. Left: 3rd order WSO 1 DIRK scheme with conventional b.c. exhibits order reduction. Middle: The same WSO 1 DIRK scheme with MBC recovers full 3rd order convergence in $u$ and $u_x$. Right: 3rd order WSO 2 DIRK scheme recovers 3rd order convergence in $u$.

**Related Publications:**- R. R. Rosales, B. Seibold, D. Shirokoff and D. Zhou.
*Order reduction in high-order implicit Runge-Kutta methods for initial boundary value problems.*Submitted to Math. Comput..

Pressure Poisson equation (PPE) reformulations of the incompressible Navier-Stokes equations
is a class of methods that replace the incompressibility constraint by a Poisson equation for the pressure, with
a proper choice of the boundary condition so that the incompressibility is maintained. This approach has important
advantages compared to the standard form:
(1) the pressure is not implicitly coupled to the velocity and can be recovered by solving a Poisson equation;
(2) no numerical boundary layers are generated, in contrast to the projection type methods.

In this work, I proposed numerical methods for solving one PPE reformulation (Shirokoff and Rosales, 2011) using mixed finite element methods and meshfree finite differences in space with high-order implicit-explicit (IMEX) time-stepping in time. This PPE reformulation is special due to its electric boundary conditions ($\hat{n}\times\vec{u}=0$ and $\nabla\cdot\vec{u}=0$). The standard nodal-based FE approach will converge to the wrong solution on domains with reentrant corners or curved boundaries. The former is due to the insufficient regularity of the solution, and the latter is a form of Babuska paradox.

The mixed finite element approach is well-suited for electric boundary conditions, and the meshfree finite difference handles curved boundaries naturally. The IMEX Runge-Kutta methods are used to treat the viscous term implicitly and the pressure explicitly. The method achieves 3rd order convergence in time with a 3rd order IMEX RK method, for problems with time-independent boundary conditions and forcing at the boundary (a special case where order reduction does not occur).

Some Benchmark Test Results (Lid-driven Cavity, Backward-facing Step)

In this work, I proposed numerical methods for solving one PPE reformulation (Shirokoff and Rosales, 2011) using mixed finite element methods and meshfree finite differences in space with high-order implicit-explicit (IMEX) time-stepping in time. This PPE reformulation is special due to its electric boundary conditions ($\hat{n}\times\vec{u}=0$ and $\nabla\cdot\vec{u}=0$). The standard nodal-based FE approach will converge to the wrong solution on domains with reentrant corners or curved boundaries. The former is due to the insufficient regularity of the solution, and the latter is a form of Babuska paradox.

The mixed finite element approach is well-suited for electric boundary conditions, and the meshfree finite difference handles curved boundaries naturally. The IMEX Runge-Kutta methods are used to treat the viscous term implicitly and the pressure explicitly. The method achieves 3rd order convergence in time with a 3rd order IMEX RK method, for problems with time-independent boundary conditions and forcing at the boundary (a special case where order reduction does not occur).

Some Benchmark Test Results (Lid-driven Cavity, Backward-facing Step)

**Related Publications:**- D. Zhou, B. Seibold, D. Shirokoff, P. Chidyagwai and R. R. Rosales.
*Meshfree finite differences for vector Poisson and pressure Poisson equations with electric boundary conditions.*Meshfree methods for Partial Differential Equations VII, Lecture Notes in Computational Science and Engineering, Vol. 100. Griebel, M. and Schweitzer, M.A. (ed.), Springer, 2015. - D. Zhou, B. Seibold, D. Shirokoff, and R. R. Rosales.
*High-order mixed finite element methods for a pressure Poisson equation reformulation of the Navier-Stokes equations with electric boundary conditions.*In preparation.

**Ph.D. Thesis:**- D. Zhou.
*High-order numerical methods for pressure Poisson equation reformulations of the incompressible Navier-Stokes equations.*Ph.D. Thesis, Temple University, 2014.

Jet schemes are based on evolving parts of the jet of the solution in time, and solving linear advection equations with high order of
accuracy by tracking characteristic curves and using suitable Hermite interpolants. Click
here for more on Jet schemes for linear advection problems.

My work is focusing on extending this scheme to nonlinear Hamilton-Jacobi type equations, where the characteristics may collide or emenate from a single point. Solutions to such problems can have kinks even with smooth initial conditions. Tracking characteristics gives rise to minimization problems at each time step. The overall scheme based on bi-cubic Hermite interpolation, approximates the entropy solution correctly and achieves optimal 3rd order accuracy in the smooth parts of the solution, while the low order approximations at the non-smooth parts remain localized.

The animation on the right is a numerical simulation of front propagation with an obstacle in the computational domain. The front is represented implicitly by a levelset function.

My work is focusing on extending this scheme to nonlinear Hamilton-Jacobi type equations, where the characteristics may collide or emenate from a single point. Solutions to such problems can have kinks even with smooth initial conditions. Tracking characteristics gives rise to minimization problems at each time step. The overall scheme based on bi-cubic Hermite interpolation, approximates the entropy solution correctly and achieves optimal 3rd order accuracy in the smooth parts of the solution, while the low order approximations at the non-smooth parts remain localized.

The animation on the right is a numerical simulation of front propagation with an obstacle in the computational domain. The front is represented implicitly by a levelset function.