Skip to content

Commit 032176e

Browse files
Merge upstream master into NonlinearSolveBase floor
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
2 parents 82521c8 + aa895b5 commit 032176e

7 files changed

Lines changed: 132 additions & 2 deletions

File tree

docs/Project.toml

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,7 @@ OrdinaryDiffEqIMEXMultistep = "9f002381-b378-40b7-97a6-27a27c83f129"
2020
OrdinaryDiffEqLinear = "521117fe-8c41-49f8-b3b6-30780b3f0fb5"
2121
OrdinaryDiffEqLowOrderRK = "1344f307-1e59-4825-a18e-ace9aa3fa4c6"
2222
OrdinaryDiffEqLowStorageRK = "b0944070-b475-4768-8dec-fb6eb410534d"
23+
OrdinaryDiffEqNewmark = "d204908a-63b9-11ef-18f5-cfec7123a93b"
2324
OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
2425
OrdinaryDiffEqNordsieck = "c9986a66-5c92-4813-8696-a7ec84c806c8"
2526
OrdinaryDiffEqPDIRK = "5dd0a6cf-3d4b-4314-aa06-06d4e299bc89"
@@ -55,6 +56,7 @@ OrdinaryDiffEqIMEXMultistep = {path = "../lib/OrdinaryDiffEqIMEXMultistep"}
5556
OrdinaryDiffEqLinear = {path = "../lib/OrdinaryDiffEqLinear"}
5657
OrdinaryDiffEqLowOrderRK = {path = "../lib/OrdinaryDiffEqLowOrderRK"}
5758
OrdinaryDiffEqLowStorageRK = {path = "../lib/OrdinaryDiffEqLowStorageRK"}
59+
OrdinaryDiffEqNewmark = {path = "../lib/OrdinaryDiffEqNewmark"}
5860
OrdinaryDiffEqNonlinearSolve = {path = "../lib/OrdinaryDiffEqNonlinearSolve"}
5961
OrdinaryDiffEqNordsieck = {path = "../lib/OrdinaryDiffEqNordsieck"}
6062
OrdinaryDiffEqPDIRK = {path = "../lib/OrdinaryDiffEqPDIRK"}
@@ -92,6 +94,7 @@ OrdinaryDiffEqIMEXMultistep = "2.0.0"
9294
OrdinaryDiffEqLinear = "2.0.0"
9395
OrdinaryDiffEqLowOrderRK = "2.0.0"
9496
OrdinaryDiffEqLowStorageRK = "3.0.0"
97+
OrdinaryDiffEqNewmark = "2.0.0"
9598
OrdinaryDiffEqNonlinearSolve = "2"
9699
OrdinaryDiffEqNordsieck = "2.0.0"
97100
OrdinaryDiffEqPDIRK = "2.0.0"

docs/make.jl

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,7 @@ using OrdinaryDiffEqNordsieck
2626
using OrdinaryDiffEqPDIRK
2727
using OrdinaryDiffEqPRK
2828
using OrdinaryDiffEqQPRK
29+
using OrdinaryDiffEqNewmark
2930
using OrdinaryDiffEqRKN
3031
using OrdinaryDiffEqRosenbrock
3132
using OrdinaryDiffEqSDIRK
@@ -69,6 +70,7 @@ makedocs(
6970
OrdinaryDiffEqPDIRK,
7071
OrdinaryDiffEqPRK,
7172
OrdinaryDiffEqQPRK,
73+
OrdinaryDiffEqNewmark,
7274
OrdinaryDiffEqRKN,
7375
OrdinaryDiffEqRosenbrock,
7476
OrdinaryDiffEqSDIRK,

docs/pages.jl

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -34,9 +34,10 @@ pages = [
3434
"imex/StabilizedIRK.md",
3535
"imex/IMEXBDF.md",
3636
],
37-
"Dynamical ODE Explicit Solvers" => [
37+
"Dynamical ODE Solvers" => [
3838
"dynamicalodeexplicit/RKN.md",
3939
"dynamicalodeexplicit/SymplecticRK.md",
40+
"dynamicalodeexplicit/Newmark.md",
4041
],
4142
"Semilinear ODE Solvers" => [
4243
"semilinear/ExponentialRK.md",
Lines changed: 96 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,96 @@
1+
```@meta
2+
CollapsedDocStrings = true
3+
```
4+
5+
# OrdinaryDiffEqNewmark
6+
7+
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.
8+
9+
## Key Properties
10+
11+
These methods provide:
12+
13+
- **Direct integration** of second-order ODEs, including mass-matrix form
14+
- **Unconditionally stable** parameter choices for structural dynamics
15+
- **Controllable high-frequency damping** via generalized-α's (`ρ∞`)
16+
- **Newmark-β and HHT-α** as special cases of generalized-α
17+
- **Adaptive time stepping** with the Zienkiewicz–Xie local truncation error estimate (Newmark-β)
18+
19+
## When to Use Newmark Methods
20+
21+
These methods are recommended for:
22+
23+
- **Structural dynamics** and finite-element/Galerkin vibration problems
24+
- **Second-order ODEs** of the form `M a = f(u, v, t)`
25+
- **Problems needing algorithmic damping** of spurious high-frequency modes
26+
- **Mass-matrix second-order systems** arising from spatial discretizations
27+
- **Cases where RKN/symplectic methods are not the right fit** (implicit structural integrators with dissipation control)
28+
29+
## Mathematical Background
30+
31+
The generalized-α method evaluates the equations of motion at interpolated states:
32+
33+
`M * aₙ₊αₘ = f(uₙ₊αf, vₙ₊αf, tₙ₊αf)`
34+
35+
with
36+
37+
- `aₙ₊αₘ = (1 - αₘ) * aₙ₊₁ + αₘ * aₙ`
38+
- `uₙ₊αf = (1 - αf) * uₙ₊₁ + αf * uₙ`
39+
- `vₙ₊αf = (1 - αf) * vₙ₊₁ + αf * vₙ`
40+
41+
and the standard Newmark updates for `uₙ₊₁` and `vₙ₊₁`. Setting `αₘ = αf = 0` recovers Newmark-β; setting `αₘ = 0` recovers HHT-α.
42+
43+
## Solver Selection Guide
44+
45+
- **`NewmarkBeta`**: Classical Newmark-β. Default `β = 1/4`, `γ = 1/2` (average acceleration, second-order when `γ = 1/2`).
46+
- **`GeneralizedAlpha`**: Preferred when controllable high-frequency damping is needed.
47+
- **`GeneralizedAlpha(; rho_inf)`**: Recommended parameterization. `ρ∞ = 1` gives no algorithmic damping (equivalent to undamped Newmark); `ρ∞ = 0` gives maximum damping. Always unconditionally stable in [0, 1].
48+
- **`GeneralizedAlpha(; alpha_hht)`**: HHT-α convenience (`α ∈ [-1/3, 0]`). Also unconditionally stable in that range.
49+
- **`GeneralizedAlpha(αm, αf, β, γ)`**: Explicit four-parameter construction. Can break unconditional stability if parameters are chosen poorly.
50+
51+
### Unconditional stability
52+
53+
**Newmark-β** is unconditionally stable when `γ ≥ 1/2` and `β ≥ (1/4) * (γ + 1/2)^2`.
54+
55+
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.
56+
57+
**Generalized-α** is unconditionally stable when `αₘ ≤ αf ≤ 1/2` and `β ≥ (1/4) * (1/2 + αf - αₘ)^2`.
58+
59+
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.
60+
61+
## Installation
62+
63+
To be able to access the solvers in `OrdinaryDiffEqNewmark`, you must first install them using the Julia package manager:
64+
65+
```julia
66+
using Pkg
67+
Pkg.add("OrdinaryDiffEqNewmark")
68+
```
69+
70+
This will only install the solvers listed at the bottom of this page.
71+
If you want to explore other solvers for your problem,
72+
you will need to install some of the other libraries listed in the navigation bar on the left.
73+
74+
## Example usage
75+
76+
```julia
77+
using OrdinaryDiffEqNewmark
78+
function f1!(dv, v, u, p, t)
79+
dv .= -u
80+
end
81+
function f2!(du, v, u, p, t)
82+
du .= v
83+
end
84+
v0 = ones(2)
85+
u0 = zeros(2)
86+
prob = DynamicalODEProblem(f1!, f2!, v0, u0, (0.0, 5.0))
87+
sol = solve(prob, NewmarkBeta(), dt = 0.1)
88+
sol_ga = solve(prob, GeneralizedAlpha(; rho_inf = 0.8), dt = 0.1)
89+
```
90+
91+
## Full list of solvers
92+
93+
```@docs
94+
NewmarkBeta
95+
GeneralizedAlpha
96+
```

lib/OrdinaryDiffEqCore/Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "OrdinaryDiffEqCore"
22
uuid = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
33
authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
4-
version = "4.7.0"
4+
version = "4.7.1"
55

66
[deps]
77
ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"

lib/OrdinaryDiffEqCore/src/integrators/integrator_interface.jl

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -511,6 +511,14 @@ function SciMLBase.reinit!(
511511
integrator.sol = sol
512512
end
513513
end
514+
# Null u0 (e.g., MTK systems with only observed/algebraic variables and no state).
515+
# `init` coerces such a `nothing` to the working state (`Float64[]`), but `prob.u0`
516+
# (the default here) and `late_binding_update_u0_p` can still hand back `nothing`.
517+
# There is no state to reinitialize, so fall back to the integrator's current `u`,
518+
# which downstream copying/saving can handle uniformly.
519+
if u0 === nothing
520+
u0 = integrator.u
521+
end
514522
if isinplace(integrator.sol.prob)
515523
recursivecopy!(integrator.u, u0)
516524
recursivecopy!(integrator.uprev, integrator.u)

test/Integrators_I/reinit_test.jl

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -147,3 +147,23 @@ end
147147
# The solver should handle the discontinuities correctly
148148
@test integrator.sol.retcode == ReturnCode.Success
149149
end
150+
151+
@testset "reinit with null u0 (no state variables)" begin
152+
# Systems with only observed/algebraic variables (e.g. MTK models that
153+
# structurally simplify to zero dynamical states) have `prob.u0 === nothing`.
154+
# `init` coerces that to an empty working state, but `reinit!` used to feed the
155+
# raw `nothing` into `recursivecopy!`/`copyat_or_push!`. A callback forces the
156+
# real ODEIntegrator path rather than the NullODEIntegrator shortcut.
157+
prob = ODEProblem((du, u, p, t) -> nothing, nothing, (0.0, 1.0))
158+
cb = DiscreteCallback((u, t, integrator) -> false, integrator -> nothing)
159+
160+
integrator = init(prob, Tsit5(); callback = cb)
161+
@test integrator.u == Float64[]
162+
solve!(integrator)
163+
@test integrator.sol.retcode == ReturnCode.Success
164+
165+
reinit!(integrator)
166+
@test integrator.u == Float64[]
167+
solve!(integrator)
168+
@test integrator.sol.retcode == ReturnCode.Success
169+
end

0 commit comments

Comments
 (0)