Coverage for pyxfem/core/nonlinear_solver.py: 49%
121 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-24 07:17 +0000
1# author Chevaugeon Nicolas 2024
2"""
3Nonlinear solvers and arc-length continuation routines.
5This module provides numerical algorithms for solving nonlinear systems of equations
6and tracing equilibrium paths in structural and continuum mechanics:
8- Newton-Raphson solvers with customizable linear solvers, finite-difference fallback
9 for Jacobians, and step-by-step convergence controls.
10- Unified root-finding interface dispatching to internal solvers or SciPy.
11- High-order arc-length (continuation) predictor routines based on Taylor expansions
12 and adaptive step-size estimation for path following.
13"""
15import copy
16from operator import isub
17from typing import TypeAlias
18from warnings import warn
20import numpy as np
21import scipy as sp
23import pyxfem.core.sparse_tools as spt
24from pyxfem.core.config import FloatArray
26InitialGuess:TypeAlias = FloatArray|tuple[FloatArray, FloatArray]|list[FloatArray]
28def _fd_jac(fun):
29 """
30 Construct a numerical Jacobian function using forward finite differences.
32 Parameters
33 ----------
34 fun : callable
35 Vector-valued function of the form ``fun(x)`` returning an array-like.
37 Returns
38 -------
39 callable
40 Function returning the 2D Jacobian matrix evaluated at ``x``.
42 Notes :
43 For 1D vector fun: ensure result is reshaped to (1, 1)
44 For N-D vector fun: approx_fprime per component or use scipy.optimize.approx_fprime
45 """
47 def jacfun(x):
48 val = sp.optimize.approx_fprime(x, fun)
49 if val.ndim == 1:
50 return val.reshape(-1, len(x))
51 return val
52 return jacfun
55def newton(fun, x0:InitialGuess, jac = None,
56 linsolver = None,
57 norm = np.linalg.norm,
58 inplace_sub = isub,
59 rtol:float=1.e-6, atol:float=1.e-6, nitermax = 10,
60 fail_on_increasing_res = False,
61 verbose = True
62 )-> sp.optimize.OptimizeResult :
63 """
64 Solve a non-linear system of equations f(x) = 0 using Newton's method.
66 Parameters
67 ----------
68 fun : callable
69 Objective function returning the residual vector at ``x``.
70 x0 : FloatArray, tuple of FloatArray, or list of FloatArray
71 Initial guess for the solution.
72 jac : callable, optional
73 Jacobian of ``fun``. If ``None``, a finite-difference approximation
74 is used.
75 linsolver : callable, optional
76 Linear solver of the form ``linsolver(A, b)``. Defaults to
77 ``numpy.linalg.solve``.
78 norm : callable, optional
79 Norm used to compute the scalar residual norm. Defaults to
80 ``numpy.linalg.norm``.
81 inplace_sub : callable, optional
82 In-place subtraction operator used to update the solution vector.
83 Defaults to ``operator.isub``.
84 rtol : float, optional
85 Relative tolerance for convergence criteria (default is 1e-6).
86 atol : float, optional
87 Absolute tolerance for convergence criteria (default is 1e-6).
88 nitermax : int, optional
89 Maximum number of Newton iterations (default is 10).
90 fail_on_increasing_res : bool, optional
91 If True, terminates iteration early if the residual increases.
92 Defaults to False.
93 verbose : bool, optional
94 Whether to print convergence progress at each iteration (default is True).
96 Returns
97 -------
98 scipy.optimize.OptimizeResult
99 Optimization result containing:
100 - ``x``: Solution array.
101 - ``success``: Boolean flag indicating successful convergence.
102 - ``message``: Termination description.
103 - ``nit``: Number of iterations performed.
104 - ``nfev``: Number of function evaluations.
105 - ``njev``: Number of Jacobian evaluations.
106 """
107 res = fun(x0)
108 nres0 = norm(res)
109 it = 0
110 nres = nres0
111 nres_prev = nres+1.
112 x = copy.deepcopy(x0)
113 if verbose:
114 print(f" {it} {nres}")
115 if jac is None :
116 warn("You should pass a jacobian function (jac) for better performances", \
117 category = UserWarning, stacklevel =2)
118 jacfun = _fd_jac(fun)
119 else:
120 jacfun = jac
121 if linsolver is None:
122 linsolver = np.linalg.solve
123 while (nres > rtol*nres0 and nres > atol and it < nitermax \
124 and ( not(fail_on_increasing_res) or (nres<nres_prev) )):
125 dx = linsolver(jacfun(x), res)
126 inplace_sub(x, dx)
127 nres_prev = nres
128 res = fun(x)
129 nres = norm (res)
130 it = it +1
131 if verbose :
132 print(f" {it} {nres}")
133 sol = sp.optimize.OptimizeResult()
134 sol.x =x
135 sol.nit = it
136 sol.nfev = it+1
137 sol.njev = it
139 if (nres<= rtol*nres0 or nres<= atol) and ( not(fail_on_increasing_res) or (nres<nres_prev)):
140 sol.success = True
141 sol.message =f'newton converged after {it:d} iterations'
142 return sol
144 sol.success = False
145 if fail_on_increasing_res and nres>=nres_prev:
146 sol.message =f'residual increased at iter {it:d} iterations'
147 return sol
148 sol.message =f'newton did not converge after nmax {nitermax:d} iterations'
149 return sol
151def root(fun, x0, args=(), method='newton', jac=None, tol=None, callback=None, options=None):
152 """
153 Find a root of a vector function.
155 Provides a unified interface dispatching either to the custom Newton solver
156 or to ``scipy.optimize.root``.
158 Parameters
159 ----------
160 fun : callable
161 Vector function to find a root of.
162 x0 : ndarray
163 Initial guess.
164 args : tuple, optional
165 Extra arguments passed to the objective function (used when method != 'newton').
166 method : str, optional
167 Solver type. Defaults to ``'newton'``. Any other value is delegated
168 to ``scipy.optimize.root``.
169 jac : bool or callable, optional
170 Jacobian evaluation function.
171 tol : float, optional
172 Tolerance for termination. Defaults to 1e-6 for ``method='newton'``.
173 callback : callable, optional
174 Callback called after each iteration (passed to SciPy solver).
175 options : dict, optional
176 Dictionary of solver options (passed to SciPy solver).
178 Returns
179 -------
180 scipy.optimize.OptimizeResult
181 The solution represented as an ``OptimizeResult`` object.
182 """
183 if method == 'newton':
184 if tol is None :
185 tol = 1.e-6
186 return newton(fun, x0, jac, rtol=tol, atol=tol, nitermax = 10)
187 else :
188 return sp.optimize.root(fun, x0, args=args, method=method, jac=jac, \
189 tol=tol, callback=callback, options=options)
192def arc_lenght_predictor(drdx, drdl, fd2rdx2_u, d2rdxdl,
193 du_prev = None, dl_prev = None, pred_order = 2,
194 alpha =None, a = None, delta = None, afac=None, amax = 0.1):
195 """
196 Compute Taylor expansion terms of the arc-length path with respect to the arc parameter.
198 Parameters
199 ----------
200 drdx : ndarray or scipy.sparse matrix
201 Tangent operator dR/dx.
202 drdl : ndarray
203 Load tangent vector dR/dlambda.
204 fd2rdx2_u : callable
205 Function computing the directional second derivative of residual
206 with respect to state: ``fd2rdx2_u(u=du)``.
207 d2rdxdl : ndarray or scipy.sparse matrix
208 Mixed second derivative d^2R/(dx dlambda).
209 du_prev : ndarray, optional
210 State increment from previous step.
211 dl_prev : float, optional
212 Load increment from previous step.
213 pred_order : {1, 2}, optional
214 Expansion order (default is 2).
215 alpha : float, optional
216 Metric parameter weighting the load factor increment. Computed automatically
217 if ``None``.
219 Returns
220 -------
221 dus : list of ndarray
222 List containing the 1st (and 2nd if ``pred_order == 2``) derivatives of state.
223 dls : list of float
224 List containing the 1st (and 2nd if ``pred_order == 2``) derivatives of load parameter.
225 alpha : float
226 Metric scaling parameter.
228 Raises
229 ------
230 ValueError
231 If ``pred_order`` is not 1 or 2.
232 """
233 dus, dls, alpha = arc_lenght_predictor_expansion(drdx, drdl, fd2rdx2_u, d2rdxdl,
234 du_prev = du_prev, dl_prev = dl_prev, pred_order = pred_order, alpha = alpha)
235 if a is None:
236 a = arc_lenght_predictor_a(dus, dls, alpha, delta, afac, amax)
237 dx, dl = arc_lenght_predictor_eval(dus, dls, a)
238 return dx, dl, alpha, a
240def arc_lenght_predictor_expansion(drdx, drdl, fd2rdx2_u, d2rdxdl,
241 du_prev = None, dl_prev = None, pred_order = 2, alpha =None):
242 """
243 Compute optimal step size 'a' for a 2nd-order arc-length expansion.
245 Parameters
246 ----------
247 du_terms : list of ndarray
248 List containing ``[du1, du2]`` state derivative vectors.
249 dl_terms : list of float
250 List containing ``[dl1, dl2]`` load parameter derivative scalars.
251 alpha : float
252 Arc-length scaling factor.
253 delta : float
254 Error tolerance parameter.
255 afac : float
256 Scale factor applied to the computed step size.
257 amax : float
258 Maximum allowable step size.
260 Returns
261 -------
262 float
263 Computed step size ``a`` bounded by ``amax``.
265 Raises
266 ------
267 ValueError
268 If ``du_terms`` does not contain exactly 2 terms.
269 """
270 # pylint: disable=invalid-name
271 K = drdx
272 q = -drdl
273 iK = spt.factorized(K)
274 ut = iK(q)
275 n = len(q)
277 if alpha is None:
278 alpha = ut.dot(ut)/n
280 def alpha_scal(u1, l1, u2, l2):
281 return u1.dot(u2)/n +alpha*l1*l2
283 def alpha_nrm2(u,l):
284 return np.sqrt(alpha_scal(u, l, u, l))
286 dl1_mag = 1./alpha_nrm2(ut,1.)
288 if du_prev is None:
289 # First step: assume forward loading
290 dl1 = dl1_mag
291 else:
292 # Subsequent steps: ensure we move forward along the path
293 # We want (dp * ut) . du_prev > 0
294 sign_dl1 = np.sign(alpha_scal(ut, dl1_mag, du_prev, dl_prev))
295 if sign_dl1 == 0:
296 sign_dl1 = 1.0 # Fallback if dot product is exactly zero (highly unlikely)
297 dl1 = sign_dl1 * dl1_mag
299 du1 = dl1*ut
300 if pred_order == 1:
301 return [du1], [dl1], alpha
302 if pred_order == 2:
303 f2_struct = -0.5*fd2rdx2_u( u = du1)
304 Kext = -d2rdxdl
305 f2 = dl1*Kext@du1 + f2_struct
306 v2 = iK(f2)
307 dl2 = -(ut.dot(v2) / n)*dl1_mag**2
308 du2 = dl2*ut +v2
309 else :
310 raise ValueError("pred_order should be 1 or 2")
311 return [du1, du2], [dl1, dl2], alpha
313def arc_lenght_predictor_a(du_terms, dl_terms, alpha, delta, afac, amax):
314 """
315 Compute optimal step size 'a' for a 2nd-order arc-length expansion.
317 Parameters
318 ----------
319 du_terms : list of ndarray
320 List containing ``[du1, du2]`` state derivative vectors.
321 dl_terms : list of float
322 List containing ``[dl1, dl2]`` load parameter derivative scalars.
323 alpha : float
324 Arc-length scaling factor.
325 delta : float
326 Error tolerance parameter.
327 afac : float
328 Scale factor applied to the computed step size.
329 amax : float
330 Maximum allowable step size.
332 Returns
333 -------
334 float
335 Computed step size ``a`` bounded by ``amax``.
337 Raises
338 ------
339 ValueError
340 If ``du_terms`` does not contain exactly 2 terms.
341 """
342 if len(du_terms) != 2:
343 raise ValueError("arc_lenght_predictor_a only work for 2nd order expansion")
344 du1, du2 = du_terms
345 dl1, dl2 = dl_terms
346 n = len(du1)
347 def alpha_scal(u1, l1, u2, l2):
348 return u1.dot(u2)/n +alpha*l1*l2
349 def alpha_nrm2(u,l):
350 return np.sqrt(alpha_scal(u, l, u, l))
351 ny1 = alpha_nrm2(du1, dl1)
352 ny2 = alpha_nrm2(du2, dl2)
353 return min((delta *ny1/ny2/ny2)**(1./3)*afac, amax)
355def arc_lenght_predictor_eval(du_terms, dl_terms, a):
356 """
357 Evaluate total state and load increments using a polynomial expansion.
359 Parameters
360 ----------
361 du_terms : list of ndarray
362 First-order (and optionally second-order) state derivative vectors.
363 dl_terms : list of float
364 First-order (and optionally second-order) load derivative scalars.
365 a : float
366 Step size parameter.
368 Returns
369 -------
370 du : ndarray
371 Combined state increment: ``a * du1 + a^2 * du2`` (if order 2).
372 dl : float
373 Combined load increment: ``a * dl1 + a^2 * dl2`` (if order 2).
374 """
375 n = len(du_terms)
376 assert 1<= n <= 2
377 du = a*du_terms[0]
378 dl = a*dl_terms[0]
379 if n == 1:
380 return du, dl
381 du+= a*a*du_terms[1]
382 dl+= a*a*dl_terms[1]
383 return du, dl