analyze stability on a nonuniform grid

Question

Assume you have a stability constraint between the space distance in time and space, for example, with an explicit Euler method for $u_t=u_{xx}$ we know $\tau\leq h^2/2$. That is, one can do stability analysis for a uniform grid, obtain a constraint and then take the minimum over a nonuniform grid for the spacing to set a constraint for the spacing in time in the case of a non-uniform grid.

Can I use the same argument for "unconditionally stable" method on non-uniform grid? My goal is to obtain a stability in the case of a nonuniform grid. Is that sufficient to consider uniform grid for the method, such as Crank-Nicolson and somehow extend the argument to the non uniform grid? What are my options obtaining stability in this case?

Answer 1

It all boils down to finding the largest or smallest eigenvalue of the stiffness matrix. You can do that analytically for a uniform mesh, and you then arrive at the traditional formula $\tau \le ch^2$. For non-uniform meshes, the relationship between mesh size and eigenvalue is not obvious any more; in particular, there is no analytic formula for the eigenvalue any more, and there is of course no longer a single mesh size $h$ either.

In other words, I'm not familiar with an approach to come up with a definite bound on the time step that makes the scheme stable in the case of non-uniform meshes. That said, the formula everyone uses is to replace $\tau \le ch^2$ by the generalization $\tau \le c \min_K h_K^2$ where $h_K$ is the mesh size of cell $K$ and you simply take the minimum over all mesh sizes. That appears to be working well for practical cases. As I said, I don't know if there's a formal proof that this works -- my suspicion is that there is no such proof for general 2d/3d meshes but I'm willing to put down a conjecture that the statement is nevertheless true.

Answer 2

Your question may have multiple answers, depending on context.

I will give a simple one, which completes and clarify the previous answer by Wolfang Bangerth, so I'm not claiming the bounty.

Within a FEM approach you may in a first step discretize only the spatial domain: \begin{equation} u^h(x,t) = \sum_i \eta_i(x) U_i(t) \end{equation} so that \begin{equation} \frac{\partial u(x,t)}{\partial t} \approx \frac{\partial u^h(x,t)}{\partial t} = \sum_i \eta_i(x) \dot U_i(t) \end{equation} where $\dot U_i$ are the time derivatives of the nodal unknowns. Applying the classical FEM method we end up, for linear problems, with an ODE's system like \begin{equation} A \dot U + K U = b \end{equation} Now we integrate with respect to time, forgetting the original PDE and the underlying mesh.

With this simple approach, stability with respect to time increment is defined at the ODE level, so that in the case of unconditionally stable methods we are not concerned about mesh. (Or at least, very bad meshes give raise to badly conditioned $K$ matrices, but this is another problem.)

In the case of explicit time integration, stability is linked to a generalized eigenvalue problem involving $A$ and $B$. Here for uniform meshes we have simple relations between $h$ and the stability limit, for non uniform meshes only estimates and bounds.

I have some experience in the field of non-linear continuum mechanics, where the ODE to be solved is \begin{equation} M(t)\ddot U + f_i(U, \dot U, \ldots \text{state variables}) = f_e(t) \end{equation} Here the stable time increment at time $t_0$ is linked to highest natural frequency $\omega_\text{max}$ (eigenvalue) of the linearized problem \begin{equation} M(t_0) \ddot U + K_{T} U = 0, \qquad U = U_\alpha \sin(\omega_\alpha t) \end{equation} where $K_T$ is the tangent stiffness matrix. Finding a good (but fast to compute) estimator of $\omega_\text{max}$ is crucial for efficiency and far from trivial.

This treatment (which is the most common in commercial non linear continuum mechanics FEM codes) does not directly addresses the problem of convergence in the time-space domain, since it relies on the simplistic assumption that convergence for the stationary problem ($t = t_0 = \text{constant}$) and a good time integrator solve the problem.

Of course the time-space convergence problem has been addressed, but as said above, this is only a partial answer to the question above

Answer 3

Let's do the analysis by discretizing in space first, then in time (method of lines). First, we take some sequence of spatial discretizations $L^i, i=\{1,2,\dotsc\}$ approximating the continuous operator $L:u \to u_{xx}$. For the sequence of discretizations $L^i$ to be a convergent discretization, we need that the sequence be consistent and stable. Since the continuum operator $L$ has all negative eigenvalues, this implies that all eigenvalues of each discrete operator $L^i$ will be in the left half plane.

Now we turn our attention to the temporal discretization of

$$u_t^i = L^i u^i$$

where all eigenvalues of $L^i$ are in the left half plane. But $A$-stable methods are defined by exactly this property. ($A$-stability is just a more precise/less misleading term for unconditional stability. There are many other kinds of stability, e.g. $L$-stability and forms of nonlinear stability.)

Note that we have not used any property of grid spacing, time step size, or even used the properties of the original PDE except that it must be stable.