# Dynamic linear elasticity problem is stiff (numerically)

I am considering a dynamic linear elasticity problem applied to a simple structure such as a beam. In system form, the PDE can be written as $$M \ddot{X} + D \dot{X} + XK = F$$ where $$X$$ represents the vector of deflections along the $$x$$ and $$y$$ directions of various points on the beam.

According to one repository I found (link here), a standard way to solve this dynamical problem is via the Newmark scheme.

Naively, why couldn't you convert this to a first-order system such as $$A \dot{Y} + BY = C$$ and solve this using backward/forward Euler?

Or is it well-known that dynamic linear elasticity problems are stiff?

There is absolutely nothing wrong with converting the second-order system to first-order form and then using appropriate numerical methods to solve it. Both implicit and explicit Euler methods can be used. However, both have only first-order accuracy, i.e. the error is proportional to the time step size. And, of course, explicit Euler is not stable unless the time step is sufficiently small.

A method with second-order accuracy can be obtained by writing the solution of

$$\frac{dy}{dt} = f(y,t)$$

as

$$y_{n+1} = y_n + h [\theta ( t_n, y_n) + (1 - \theta ) f ( t_{n+1}, y_{n+1})]$$ where $$h$$ is the time step, $$h=t_{n+1}-t_n$$. Selecting $$\theta=0$$ gives implicit Euler, $$\theta=1$$ gives explicit Euler, and $$\theta=1/2$$ gives the second-order accurate trapezoidal method.

The Newmark numerical "method" for solving second-order ODE is actually a family of methods where a specific method can be chosen by selecting values for the $$\beta$$ and $$\gamma$$ parameters. One of the more popular methods in this family is obtained by selecting $$\beta=1/4$$ and $$\gamma=1/2$$ which gives a constant acceleration and a linear velocity over the time step. Using this method on the second-order equations gives identical results to using the trapezoidal method on the first order version of the equations.

I would not say the that elasticity equations are inherently stiff. The stiffness results when they are discretized in space using FEM or some other approach. The system of ODE becomes stiffer as the spatial mesh is refined. This is why implicit ODE methods are often used for these types of problems.

• Actually, the original elasticity equations can be said to be infinitely stiff because the ratio of the largest to the smallest eigenvalue of the elastic operator is infinite. – Wolfgang Bangerth Nov 26 '19 at 17:13
• What does it even mean to talk about "stiffness" of the underlying PDE? I guess I should have been specific about that in my answer. – Bill Greene Nov 26 '19 at 18:14
• "Stiff" is commonly interpreted as a system having a large ratio of shortest to longest time scales, which corresponds to the ratio of largest to smallest eigenvalues of the Jacobian of the right hand side of an ODE. Here, that right hand side corresponds to the matrix $K$ (plus lower order terms that do not affect the ratio). The smallest eigenvalue is going to be $O(1)$, the largest of order $O(h^{-2})$ (or $O(h^2)$ and $O(1)$), yielding a degree of ill-posedness of $O(h^{-2})$. My statement above relates to the eigenvalues of the continuous Laplace operator. – Wolfgang Bangerth Nov 26 '19 at 22:13