# analyze stability on a nonuniform grid

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?

• I think I don't quite understand what you're looking for. For unconditionally stable methods, there is no stability condition such as $\tau \le h^2/2$ -- unconditionally stable means that there is no condition. Or are you asking how to derive stability estimates for non-uniform meshes? – Wolfgang Bangerth Jul 16 '12 at 18:43
• To claim the the method is unconditionally stable means I start with an assumption that the grid is uniform and do the analysis. However, if the grid in non uniform I can't use the same proof and to do the proof in a latter case seems to be much harder! – Kamil Jul 16 '12 at 22:22
• But then the term "unconditional" doesn't make any sense. "Unconditional" means that there is no condition that has to be satisfied for the scheme to be stable. – Wolfgang Bangerth Jul 16 '12 at 23:23

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.

• well, the fact is that there are cases when people use nonuniform grid. I can do the error estimates and see that in the case of a uniform grid the numerical results confirm theoretical results. However, I can see that I have the same order of convergence when I use nonuniform grid, the question is would I need to prove that the method retains the same accuracy once I go from a uniform grid to nonuniform or there is some statement that says it is enough to show error estimates for a uniform and then.... – Kamil Jul 17 '12 at 0:54
• No, I don't think you can just generalize to non-uniform meshes. I mean, estimating the spatial error on nonuniform meshes is not substantially more difficult than on uniform meshes at least for the finite element method. But proving that a time stepping scheme is stable for a particular size of the time step $\tau$ may not be so trivial since it involves finding the largest eigenvalue of the stiffness matrix and as I mentioned I'm not aware of any analytical formulas for that. – Wolfgang Bangerth Jul 17 '12 at 9:20
• Well, the proof for stability by energy methods for finite difference method such as Crank-Nicolson is based on summation by parts and therefore it gets harder as we go from uniform to nonuniform grid as, for example, $(\partial^+u,v)=-(u,\partial^-v)$ doesn't hold anymore, where the above are forward and backward difference operator. – Kamil Jul 17 '12 at 13:53
• Oh, I might have misunderstood your question from the beginning. Are you talking about a non-uniform spatial grid, or a non-uniform temporal grid? – Wolfgang Bangerth Jul 17 '12 at 20:01
• non uniform in space, i.e. variable $x$ for $u_t=u_{xx}$ – Kamil Jul 18 '12 at 0:55

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: $$u^h(x,t) = \sum_i \eta_i(x) U_i(t)$$ so that $$\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)$$ 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 $$A \dot U + K U = b$$ 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 $$M(t)\ddot U + f_i(U, \dot U, \ldots \text{state variables}) = f_e(t)$$ Here the stable time increment at time $t_0$ is linked to highest natural frequency $\omega_\text{max}$ (eigenvalue) of the linearized problem $$M(t_0) \ddot U + K_{T} U = 0, \qquad U = U_\alpha \sin(\omega_\alpha t)$$ 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

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.

• So, are you saying that for any type of grid I have eigenvalues of $L_i$ with non-positive real part? So, as long as the mesh is not crazy and the approximation for the stable operator $L$ is consistent I have $A$-steble method? I guess eigenvalues for $L_i$ change with different mesh but stay on the left side? – Kamil Jul 17 '12 at 13:49
• @Jed, there is one problem (at least in theory) with your argument. Depending on the boundary conditions, the Laplacian may have a zero eigenvalue. A consistent spatial discretization could therefore have an eigenvalue with positive real part and exhibit a sort of instability. It would still be Lax-Richtmeyer stable, though, in an appropriate sense. – David Ketcheson Jul 17 '12 at 22:22
• Also, it should be noted that this argument would not apply in the case of a hyperbolic PDE. – David Ketcheson Jul 17 '12 at 22:23