Transient analysis
==================

:func:`femorph_solver.solvers.transient.solve_transient` integrates the
structural equation of motion

.. math::

    M \, \ddot{u} + C \, \dot{u} + K \, u = F(t)

over a uniform timestep grid using the Newmark-β scheme.  With the
default parameters :math:`\beta = 1/4,\ \gamma = 1/2` the integrator
reduces to the unconditionally-stable average-acceleration method; any
consistent pair ``(β, γ)`` can be used.

Algorithm
---------

Newmark-β updates:

.. math::

    u_{n+1} &= u_n + \Delta t\, \dot{u}_n
              + \Delta t^{2}\bigl[\tfrac{1}{2}(1-2\beta)\ddot{u}_n
              + \beta\,\ddot{u}_{n+1}\bigr] \\
    \dot{u}_{n+1} &= \dot{u}_n + \Delta t\bigl[(1-\gamma)\ddot{u}_n
              + \gamma\,\ddot{u}_{n+1}\bigr]

Rearranging the equation of motion at :math:`t_{n+1}` gives an
effective linear system

.. math::

    A \, u_{n+1} = b_{n+1}

with :math:`A = a_0 M + a_1 C + K` and
:math:`a_0 = 1/(\beta \Delta t^2)`,
:math:`a_1 = \gamma / (\beta \Delta t)`.
:math:`A` is factored **once** outside the loop; every subsequent
time step is a cheap back-solve.

API
---

.. code-block:: python

    from femorph_solver.solvers.transient import solve_transient

    K = m.stiffness_matrix()
    M = m.mass_matrix()

    # Constant step load (no damping):
    result = solve_transient(
        K, M, F=load_vector,
        dt=1e-4,
        n_steps=1000,
        prescribed={d: 0.0 for d in fixed_dofs},
    )

    result.time             # (n_steps + 1,) time grid
    result.displacement     # (n_steps + 1, N)
    result.velocity         # (n_steps + 1, N)
    result.acceleration     # (n_steps + 1, N)

Time-varying loads pass a callable instead of an array::

    def load(t):
        return force_template * np.sin(2 * np.pi * 50 * t)

    result = solve_transient(K, M, F=load, dt=1e-4, n_steps=2000,
                             prescribed=prescribed)

Damping
-------

Pass a sparse damping matrix ``C`` if you have one::

    C = 0.01 * M + 0.001 * K       # Rayleigh damping
    result = solve_transient(K, M, F=load, dt=1e-4, n_steps=1000, C=C,
                             prescribed=prescribed)

Without a ``C``, the integrator runs the undamped form.  femorph-solver
does not (yet) synthesise damping from material records; build the
matrix yourself with whatever model is appropriate.

Initial conditions
------------------

``u0`` and ``v0`` set the starting displacement / velocity.  Omitting
them zero-initialises both.  The initial acceleration is derived from
the equation of motion at :math:`t=0` so ``M ü₀ = F(0) - K u₀ - C v₀``
is satisfied exactly.

Explicit central-difference
---------------------------

Setting ``beta=0, gamma=0.5`` would give the explicit central-difference
method, but the current implementation factors :math:`A` with
``scipy.sparse.linalg.splu`` and thus requires :math:`\beta > 0`.
Explicit paths need a separate mass-factoring code path not yet
implemented.

Thread control
--------------

``thread_limit`` wraps both the one-off :math:`A` factorisation and
every back-solve.  See :doc:`threads`.

Stability notes
---------------

- **Unconditional stability** only applies to :math:`\gamma \geq 1/2`
  and :math:`\beta \geq (1/4)(1/2 + \gamma)^2`.  The default
  :math:`(1/4, 1/2)` sits exactly on that bound.
- **Numerical damping** appears if :math:`\gamma > 1/2`.  The
  trapezoid variant :math:`(1/4, 1/2)` is energy-conservative on
  linear problems.

See also
--------

- :doc:`static` — the steady-state limit :math:`(\dot u = \ddot u = 0)`.
- :doc:`modal` — free-vibration eigenmodes of :math:`(K, M)`.
- :doc:`choosing-a-solver` — which sparse backend factors :math:`A`.