# author Chevaugeon Nicolas 2024
"""
Nonlinear solvers and arc-length continuation routines.
This module provides numerical algorithms for solving nonlinear systems of equations
and tracing equilibrium paths in structural and continuum mechanics:
- Newton-Raphson solvers with customizable linear solvers, finite-difference fallback
for Jacobians, and step-by-step convergence controls.
- Unified root-finding interface dispatching to internal solvers or SciPy.
- High-order arc-length (continuation) predictor routines based on Taylor expansions
and adaptive step-size estimation for path following.
"""
import copy
from operator import isub
from typing import TypeAlias
from warnings import warn
import numpy as np
import scipy as sp
import pyxfem.core.sparse_tools as spt
from pyxfem.core.config import FloatArray
InitialGuess:TypeAlias = FloatArray|tuple[FloatArray, FloatArray]|list[FloatArray]
def _fd_jac(fun):
"""
Construct a numerical Jacobian function using forward finite differences.
Parameters
----------
fun : callable
Vector-valued function of the form ``fun(x)`` returning an array-like.
Returns
-------
callable
Function returning the 2D Jacobian matrix evaluated at ``x``.
Notes :
For 1D vector fun: ensure result is reshaped to (1, 1)
For N-D vector fun: approx_fprime per component or use scipy.optimize.approx_fprime
"""
def jacfun(x):
val = sp.optimize.approx_fprime(x, fun)
if val.ndim == 1:
return val.reshape(-1, len(x))
return val
return jacfun
[docs]
def newton(fun, x0:InitialGuess, jac = None,
linsolver = None,
norm = np.linalg.norm,
inplace_sub = isub,
rtol:float=1.e-6, atol:float=1.e-6, nitermax = 10,
fail_on_increasing_res = False,
verbose = True
)-> sp.optimize.OptimizeResult :
"""
Solve a non-linear system of equations f(x) = 0 using Newton's method.
Parameters
----------
fun : callable
Objective function returning the residual vector at ``x``.
x0 : FloatArray, tuple of FloatArray, or list of FloatArray
Initial guess for the solution.
jac : callable, optional
Jacobian of ``fun``. If ``None``, a finite-difference approximation
is used.
linsolver : callable, optional
Linear solver of the form ``linsolver(A, b)``. Defaults to
``numpy.linalg.solve``.
norm : callable, optional
Norm used to compute the scalar residual norm. Defaults to
``numpy.linalg.norm``.
inplace_sub : callable, optional
In-place subtraction operator used to update the solution vector.
Defaults to ``operator.isub``.
rtol : float, optional
Relative tolerance for convergence criteria (default is 1e-6).
atol : float, optional
Absolute tolerance for convergence criteria (default is 1e-6).
nitermax : int, optional
Maximum number of Newton iterations (default is 10).
fail_on_increasing_res : bool, optional
If True, terminates iteration early if the residual increases.
Defaults to False.
verbose : bool, optional
Whether to print convergence progress at each iteration (default is True).
Returns
-------
scipy.optimize.OptimizeResult
Optimization result containing:
- ``x``: Solution array.
- ``success``: Boolean flag indicating successful convergence.
- ``message``: Termination description.
- ``nit``: Number of iterations performed.
- ``nfev``: Number of function evaluations.
- ``njev``: Number of Jacobian evaluations.
"""
res = fun(x0)
nres0 = norm(res)
it = 0
nres = nres0
nres_prev = nres+1.
x = copy.deepcopy(x0)
if verbose:
print(f" {it} {nres}")
if jac is None :
warn("You should pass a jacobian function (jac) for better performances", \
category = UserWarning, stacklevel =2)
jacfun = _fd_jac(fun)
else:
jacfun = jac
if linsolver is None:
linsolver = np.linalg.solve
while (nres > rtol*nres0 and nres > atol and it < nitermax \
and ( not(fail_on_increasing_res) or (nres<nres_prev) )):
dx = linsolver(jacfun(x), res)
inplace_sub(x, dx)
nres_prev = nres
res = fun(x)
nres = norm (res)
it = it +1
if verbose :
print(f" {it} {nres}")
sol = sp.optimize.OptimizeResult()
sol.x =x
sol.nit = it
sol.nfev = it+1
sol.njev = it
if (nres<= rtol*nres0 or nres<= atol) and ( not(fail_on_increasing_res) or (nres<nres_prev)):
sol.success = True
sol.message =f'newton converged after {it:d} iterations'
return sol
sol.success = False
if fail_on_increasing_res and nres>=nres_prev:
sol.message =f'residual increased at iter {it:d} iterations'
return sol
sol.message =f'newton did not converge after nmax {nitermax:d} iterations'
return sol
[docs]
def root(fun, x0, args=(), method='newton', jac=None, tol=None, callback=None, options=None):
"""
Find a root of a vector function.
Provides a unified interface dispatching either to the custom Newton solver
or to ``scipy.optimize.root``.
Parameters
----------
fun : callable
Vector function to find a root of.
x0 : ndarray
Initial guess.
args : tuple, optional
Extra arguments passed to the objective function (used when method != 'newton').
method : str, optional
Solver type. Defaults to ``'newton'``. Any other value is delegated
to ``scipy.optimize.root``.
jac : bool or callable, optional
Jacobian evaluation function.
tol : float, optional
Tolerance for termination. Defaults to 1e-6 for ``method='newton'``.
callback : callable, optional
Callback called after each iteration (passed to SciPy solver).
options : dict, optional
Dictionary of solver options (passed to SciPy solver).
Returns
-------
scipy.optimize.OptimizeResult
The solution represented as an ``OptimizeResult`` object.
"""
if method == 'newton':
if tol is None :
tol = 1.e-6
return newton(fun, x0, jac, rtol=tol, atol=tol, nitermax = 10)
else :
return sp.optimize.root(fun, x0, args=args, method=method, jac=jac, \
tol=tol, callback=callback, options=options)
[docs]
def arc_lenght_predictor(drdx, drdl, fd2rdx2_u, d2rdxdl,
du_prev = None, dl_prev = None, pred_order = 2,
alpha =None, a = None, delta = None, afac=None, amax = 0.1):
"""
Compute Taylor expansion terms of the arc-length path with respect to the arc parameter.
Parameters
----------
drdx : ndarray or scipy.sparse matrix
Tangent operator dR/dx.
drdl : ndarray
Load tangent vector dR/dlambda.
fd2rdx2_u : callable
Function computing the directional second derivative of residual
with respect to state: ``fd2rdx2_u(u=du)``.
d2rdxdl : ndarray or scipy.sparse matrix
Mixed second derivative d^2R/(dx dlambda).
du_prev : ndarray, optional
State increment from previous step.
dl_prev : float, optional
Load increment from previous step.
pred_order : {1, 2}, optional
Expansion order (default is 2).
alpha : float, optional
Metric parameter weighting the load factor increment. Computed automatically
if ``None``.
Returns
-------
dus : list of ndarray
List containing the 1st (and 2nd if ``pred_order == 2``) derivatives of state.
dls : list of float
List containing the 1st (and 2nd if ``pred_order == 2``) derivatives of load parameter.
alpha : float
Metric scaling parameter.
Raises
------
ValueError
If ``pred_order`` is not 1 or 2.
"""
dus, dls, alpha = arc_lenght_predictor_expansion(drdx, drdl, fd2rdx2_u, d2rdxdl,
du_prev = du_prev, dl_prev = dl_prev, pred_order = pred_order, alpha = alpha)
if a is None:
a = arc_lenght_predictor_a(dus, dls, alpha, delta, afac, amax)
dx, dl = arc_lenght_predictor_eval(dus, dls, a)
return dx, dl, alpha, a
[docs]
def arc_lenght_predictor_expansion(drdx, drdl, fd2rdx2_u, d2rdxdl,
du_prev = None, dl_prev = None, pred_order = 2, alpha =None):
"""
Compute optimal step size 'a' for a 2nd-order arc-length expansion.
Parameters
----------
du_terms : list of ndarray
List containing ``[du1, du2]`` state derivative vectors.
dl_terms : list of float
List containing ``[dl1, dl2]`` load parameter derivative scalars.
alpha : float
Arc-length scaling factor.
delta : float
Error tolerance parameter.
afac : float
Scale factor applied to the computed step size.
amax : float
Maximum allowable step size.
Returns
-------
float
Computed step size ``a`` bounded by ``amax``.
Raises
------
ValueError
If ``du_terms`` does not contain exactly 2 terms.
"""
# pylint: disable=invalid-name
K = drdx
q = -drdl
iK = spt.factorized(K)
ut = iK(q)
n = len(q)
if alpha is None:
alpha = ut.dot(ut)/n
def alpha_scal(u1, l1, u2, l2):
return u1.dot(u2)/n +alpha*l1*l2
def alpha_nrm2(u,l):
return np.sqrt(alpha_scal(u, l, u, l))
dl1_mag = 1./alpha_nrm2(ut,1.)
if du_prev is None:
# First step: assume forward loading
dl1 = dl1_mag
else:
# Subsequent steps: ensure we move forward along the path
# We want (dp * ut) . du_prev > 0
sign_dl1 = np.sign(alpha_scal(ut, dl1_mag, du_prev, dl_prev))
if sign_dl1 == 0:
sign_dl1 = 1.0 # Fallback if dot product is exactly zero (highly unlikely)
dl1 = sign_dl1 * dl1_mag
du1 = dl1*ut
if pred_order == 1:
return [du1], [dl1], alpha
if pred_order == 2:
f2_struct = -0.5*fd2rdx2_u( u = du1)
Kext = -d2rdxdl
f2 = dl1*Kext@du1 + f2_struct
v2 = iK(f2)
dl2 = -(ut.dot(v2) / n)*dl1_mag**2
du2 = dl2*ut +v2
else :
raise ValueError("pred_order should be 1 or 2")
return [du1, du2], [dl1, dl2], alpha
[docs]
def arc_lenght_predictor_a(du_terms, dl_terms, alpha, delta, afac, amax):
"""
Compute optimal step size 'a' for a 2nd-order arc-length expansion.
Parameters
----------
du_terms : list of ndarray
List containing ``[du1, du2]`` state derivative vectors.
dl_terms : list of float
List containing ``[dl1, dl2]`` load parameter derivative scalars.
alpha : float
Arc-length scaling factor.
delta : float
Error tolerance parameter.
afac : float
Scale factor applied to the computed step size.
amax : float
Maximum allowable step size.
Returns
-------
float
Computed step size ``a`` bounded by ``amax``.
Raises
------
ValueError
If ``du_terms`` does not contain exactly 2 terms.
"""
if len(du_terms) != 2:
raise ValueError("arc_lenght_predictor_a only work for 2nd order expansion")
du1, du2 = du_terms
dl1, dl2 = dl_terms
n = len(du1)
def alpha_scal(u1, l1, u2, l2):
return u1.dot(u2)/n +alpha*l1*l2
def alpha_nrm2(u,l):
return np.sqrt(alpha_scal(u, l, u, l))
ny1 = alpha_nrm2(du1, dl1)
ny2 = alpha_nrm2(du2, dl2)
return min((delta *ny1/ny2/ny2)**(1./3)*afac, amax)
[docs]
def arc_lenght_predictor_eval(du_terms, dl_terms, a):
"""
Evaluate total state and load increments using a polynomial expansion.
Parameters
----------
du_terms : list of ndarray
First-order (and optionally second-order) state derivative vectors.
dl_terms : list of float
First-order (and optionally second-order) load derivative scalars.
a : float
Step size parameter.
Returns
-------
du : ndarray
Combined state increment: ``a * du1 + a^2 * du2`` (if order 2).
dl : float
Combined load increment: ``a * dl1 + a^2 * dl2`` (if order 2).
"""
n = len(du_terms)
assert 1<= n <= 2
du = a*du_terms[0]
dl = a*dl_terms[0]
if n == 1:
return du, dl
du+= a*a*du_terms[1]
dl+= a*a*dl_terms[1]
return du, dl