Regular Newmark-β Method

The update procedure of the Newmark-β method is:


where the accelerations are given by solving Mu¨+Cu˙+f(u)=0. For a single-dof toy problem where f(u)=ku and there is no damping, then we just have u¨=kmu=ω2u. Plugging this back into the update procedure above gives


Rearrange things so that all the new stuff is on the left, and the old stuff on the right


and put it in matrix form so we can actually see what's going on



If the update matrix A had an eigenvalue with magnitude greater than 1, then applying it over and over would make u,u˙ grow without bound-- which is obviously wrong. So, which combinations of {β,γ,Δt} ensure that the eigenvalues of A stay inside the unit disk?

Since A is only 2x2, we can just write out its eigenvalues in closed form (with τ:=ωΔt):



In theory, we could visualize the 3D region defined by ||λ1,2(β,γ,τ)||1, but instead we'll just fix γ=12 (since there is no point using Newmark-β with any other value of γ) and plot the stability region in 2D:





This shows the classical stability results:



Constrained Newmark-β Method

What happens if we're using Newmark-β to integrate through time, but we would also like to prescribe the motion of a degree of freedom to match some function of time, U(t)?


Displacement Control

One natural way to constrain the motion of that degree of freedom is to modify the way the acceleration is calculated to ensure that the displacement takes on the right value at the end of the step:


This, together with the velocity update step, give us two equations relating {u¨i+1,u˙i+1} to {u¨i,u˙i}, so we can figure out the update matrix for this case as well:


Like before, calculate the eigenvalues to figure out the stability criteria


Note: this eigenvalues of this update matrix do not depend on Δt, so any instabilities associated with this enforcement approach do not go away with timestep refinement!

and plot the stability region



This region is equivalent to the inequality


which is a well-known sufficient condition for Newmark-β stability, but there is an important difference in this context: it's not a sufficient condition any more, it's a necessary condition. This means that some of the common choices of β=16 (linear acceleration) and β=0 (central difference) do not work with displacement control, for any timestep.

For example, consider what happens when trying to directly prescribe U(t)=14sin(2πt), with β=16.


While it's not surprising to see a time integrator go unstable with large timesteps, it is surprising to see that taking arbitrarily small timesteps doesn't help. This goes back to the note earlier about how the stability of the constrained dofs doesn't depend on the timestep. So, displacement control cannot be used with linear acceleration or central difference methods.

Velocity Control

Another possible interpretation of the constraint is u˙i+1=U˙(ti+1). If we put this assumption into the update equations like before and turn the algebraic crank, we eventually get


This time around, the eigenvalues of the update matrix are simple


and the expression for λ just tells us that γ12 (which is already satisfied by any of the useful Newmark methods).



Mathematica notebooks are available for the stability analyses and a basic numerical example.