GitHub

OrdinaryDiffEqNewmark

Newmark-β and generalized-α methods for second-order ODEs, typically in mass-matrix form. These are second-order time integrators, advancing the displacement and velocity with Newmark updates and solving a nonlinear residual for the acceleration directly. They arose from time stepping in structural and computational mechanics and are designed for such systems.

Key Properties

These methods provide:

  • Direct integration of second-order ODEs, including mass-matrix form
  • Unconditionally stable parameter choices for structural dynamics
  • Controllable high-frequency damping via generalized-α's (ρ∞)
  • Newmark-β and HHT-α as special cases of generalized-α
  • Adaptive time stepping with the Zienkiewicz–Xie local truncation error estimate (Newmark-β)

When to Use Newmark Methods

These methods are recommended for:

  • Structural dynamics and finite-element/Galerkin vibration problems
  • Second-order ODEs of the form M a = f(u, v, t)
  • Problems needing algorithmic damping of spurious high-frequency modes
  • Mass-matrix second-order systems arising from spatial discretizations
  • Cases where RKN/symplectic methods are not the right fit (implicit structural integrators with dissipation control)

Mathematical Background

The generalized-α method evaluates the equations of motion at interpolated states:

M * aₙ₊αₘ = f(uₙ₊αf, vₙ₊αf, tₙ₊αf)

with

  • aₙ₊αₘ = (1 - αₘ) * aₙ₊₁ + αₘ * aₙ
  • uₙ₊αf = (1 - αf) * uₙ₊₁ + αf * uₙ
  • vₙ₊αf = (1 - αf) * vₙ₊₁ + αf * vₙ

and the standard Newmark updates for uₙ₊₁ and vₙ₊₁. Setting αₘ = αf = 0 recovers Newmark-β; setting αₘ = 0 recovers HHT-α.

Solver Selection Guide

  • NewmarkBeta: Classical Newmark-β. Default β = 1/4, γ = 1/2 (average acceleration, second-order when γ = 1/2).
  • GeneralizedAlpha: Preferred when controllable high-frequency damping is needed.
    • GeneralizedAlpha(; rho_inf): Recommended parameterization. ρ∞ = 1 gives no algorithmic damping (equivalent to undamped Newmark); ρ∞ = 0 gives maximum damping. Always unconditionally stable in [0, 1].
    • GeneralizedAlpha(; alpha_hht): HHT-α convenience (α ∈ [-1/3, 0]). Also unconditionally stable in that range.
    • GeneralizedAlpha(αm, αf, β, γ): Explicit four-parameter construction. Can break unconditional stability if parameters are chosen poorly.

Unconditional stability

Newmark-β is unconditionally stable when γ ≥ 1/2 and β ≥ (1/4) * (γ + 1/2)^2.

The default β = 1/4, γ = 1/2 (average acceleration) sits on this bound and is second-order. γ = 1/2 with β < 1/4 (e.g. central difference β = 0) is only conditionally stable. γ > 1/2 adds numerical damping but drops the method to first order.

Generalized-α is unconditionally stable when αₘ ≤ αf ≤ 1/2 and β ≥ (1/4) * (1/2 + αf - αₘ)^2.

Second-order accuracy further requires γ = 1/2 - αₘ + αf. The rho_inf and alpha_hht constructors enforce these choices; the four-parameter form asserts the stability inequalities above at construction, so values that violate them will error rather than silently run unstably.

Installation

To be able to access the solvers in OrdinaryDiffEqNewmark, you must first install them using the Julia package manager:

using Pkg
Pkg.add("OrdinaryDiffEqNewmark")

This will only install the solvers listed at the bottom of this page. If you want to explore other solvers for your problem, you will need to install some of the other libraries listed in the navigation bar on the left.

Example usage

using OrdinaryDiffEqNewmark
function f1!(dv, v, u, p, t)
    return dv .= -u
end
function f2!(du, v, u, p, t)
    return du .= v
end
v0 = ones(2)
u0 = zeros(2)
prob = DynamicalODEProblem(f1!, f2!, v0, u0, (0.0, 5.0))
sol = solve(prob, NewmarkBeta(), dt = 0.1)
sol_ga = solve(prob, GeneralizedAlpha(; rho_inf = 0.8), dt = 0.1)

Full list of solvers

OrdinaryDiffEqNewmark.NewmarkBeta — Type
NewmarkBeta

Classical Newmark-β method to solve second order ODEs, possibly in mass matrix form. Local truncation errors are estimated with the estimate of Zienkiewicz and Xie.

Adaptive time stepping

The Zienkiewicz-Xie estimate monitors both solution components. The position channel uses the leading local error (β - 1/6)·dt²·(aₙ₊₁ - aₙ). The velocity channel uses (γ - 1/2)·dt·(aₙ₊₁ - aₙ) with the coefficient magnitude clamped at 1/100 from below: on and near the second-order family γ = 1/2 the true coefficient vanishes, and the channel falls back to the surrogate dt·(aₙ₊₁ - aₙ)/100, one order below the true velocity error there. This keeps the velocity error tracking the tolerance at all frequencies — measured on u'' = -ω²u for ω = 1..1000, tol = 1e-4..1e-10, the global error stays below 2.5·nsteps·tol everywhere, where without the clamp the velocity error reached the full peak-to-peak range with a Success return at ω = 1000 — at the price of dt ~ tol^(1/2) instead of tol^(1/3) where the surrogate binds: step counts at ω = 1 grow by 1.05x (tol = 1e-4) up to 9.1x (tol = 1e-10). Combinations with β within 1/24 of 1/6 and γ within 1/100 of 1/2 are rejected in adaptive mode because neither leading error coefficient is measurable there; the linear-acceleration method β = 1/6, γ = 1/2 remains available with adaptive = false.

References

Newmark, Nathan (1959), "A method of computation for structural dynamics", Journal of the Engineering Mechanics Division, 85 (EM3) (3): 67–94, doi: https://doi.org/10.1061/JMCEA3.0000098

Zienkiewicz, O. C., and Y. M. Xie. "A simple error estimator and adaptive time stepping procedure for dynamic analysis." Earthquake engineering & structural dynamics 20.9 (1991): 871-887, doi: https://doi.org/10.1002/eqe.4290200907

source
OrdinaryDiffEqNewmark.GeneralizedAlpha — Type
GeneralizedAlpha

Generalized-α method for second-order ODEs in mass-matrix form, due to Chung & Hulbert (1993). Encompasses Newmark-β (αₘ = αf = 0) and HHT-α (αₘ = 0) as special cases.

The method evaluates the equations of motion at interpolated states:

M · aₙ₊αₘ = f(uₙ₊αf, vₙ₊αf, tₙ₊αf)

where:

aₙ₊αₘ = (1 - αₘ) · aₙ₊₁ + αₘ · aₙ
uₙ₊αf = (1 - αf) · uₙ₊₁ + αf · uₙ
vₙ₊αf = (1 - αf) · vₙ₊₁ + αf · vₙ

with the standard Newmark update formulas for uₙ₊₁ and vₙ₊₁.

Constructors

GeneralizedAlpha(; rho_inf)                  # spectral radius ρ∞ ∈ [0, 1]
GeneralizedAlpha(αm, αf, β, γ)               # all four parameters directly
GeneralizedAlpha(; alpha_hht)                 # HHT-α convenience (αₘ = 0)

ρ∞ parameterization (recommended)

ρ∞ ∈ [0, 1] is the spectral radius at infinity (high-frequency damping). ρ∞ = 1 → no algorithmic damping (identical to Newmark with γ = 1/2, β = 1/4). ρ∞ = 0 → maximum algorithmic damping.

Parameters are set optimally (Chung & Hulbert 1993):

αₘ = (2ρ∞ - 1) / (ρ∞ + 1)
αf  = ρ∞        / (ρ∞ + 1)
γ   = 1/2 - αₘ + αf
β   = (1/2 + αf - αₘ)² / 4

HHT-α convenience constructor

GeneralizedAlpha(; alpha_hht = -0.1)   # α ∈ [-1/3, 0]

Sets αₘ = 0, αf = -α, γ = (1 - 2α)/2, β = (1 - α)²/4.

Order and adaptive time stepping

This implementation re-derives the acceleration from f at the accepted state on every step instead of carrying the algorithmic acceleration, so it is second order iff γ(1 - αf) = (1 - αm)/2, which for αₘ ≠ αf differs from the textbook condition γ = 1/2 - αₘ + αf; the ρ∞ and HHT parameterizations with damping therefore converge at first order with a small leading constant (measured fixed-step order 2.00 on the implementation locus).

The Zienkiewicz-Xie estimate is extended with the coefficients β - 1/(6κ) (position) and γ - 1/(2κ) (velocity), κ = (1 - αf)/(1 - αm), the velocity coefficient magnitude clamped at 1/100 from below so that the velocity channel stays live on the family γ = 1/(2κ) where it otherwise vanishes. Large ρ∞ sits inside that clamp (ρ∞ = 0.9 has γ - 1/(2κ) ≈ 0.0026, ρ∞ = 1 is exactly on the family); measured on u'' = -ω²u at ω = 100, ρ∞ = 0.9 holds the global error below 0.7·nsteps·tol for tol = 1e-4..1e-6. Parameter sets with β within 1/24 of 1/(6κ) and γ within 1/100 of 1/(2κ) are rejected in adaptive mode because neither leading error coefficient is measurable there; they remain available with adaptive = false.

References

Chung, J., and Hulbert, G. M. (1993), "A time integration algorithm for structural dynamics with improved numerical dissipation: The generalized-α method", Journal of Applied Mechanics, 60(2): 371-375. doi: https://doi.org/10.1115/1.2900803

Hilber, H. M., Hughes, T. J. R., and Taylor, R. L. (1977), "Improved numerical dissipation for time integration algorithms in structural dynamics", Earthquake Engineering & Structural Dynamics, 5(3): 283-292. doi: https://doi.org/10.1002/eqe.4290050306

source