Source code for pymor.algorithms.newton

# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2019 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)

import numpy as np

from pymor.core.defaults import defaults
from pymor.core.exceptions import InversionError, NewtonError
from pymor.core.logger import getLogger


[docs]@defaults('miniter', 'maxiter', 'rtol', 'atol', 'stagnation_window', 'stagnation_threshold') def newton(operator, rhs, initial_guess=None, mu=None, error_norm=None, least_squares=False, miniter=0, maxiter=100, rtol=-1., atol=-1., stagnation_window=3, stagnation_threshold=0.9, return_stages=False, return_residuals=False): """Basic Newton algorithm. This method solves the nonlinear equation :: A(U, mu) = V for `U` using the Newton method. Parameters ---------- operator The |Operator| `A`. `A` must implement the :meth:`~pymor.operators.interfaces.OperatorInterface.jacobian` interface method. rhs |VectorArray| of length 1 containing the vector `V`. initial_guess If not `None`, a |VectorArray| of length 1 containing an initial guess for the solution `U`. mu The |Parameter| for which to solve the equation. error_norm The norm with which the norm of the residual is computed. If `None`, the Euclidean norm is used. least_squares If `True`, use a least squares linear solver (e.g. for residual minimization). miniter Minimum amount of iterations to perform. maxiter Fail if the iteration count reaches this value without converging. rtol Finish when the residual norm has been reduced by this factor relative to the norm of the initial residual. atol Finish when the residual norm is below this threshold. stagnation_window Finish when the residual norm has not been reduced by a factor of `stagnation_threshold` during the last `stagnation_window` iterations. stagnation_threshold See `stagnation_window`. return_stages If `True`, return a |VectorArray| of the intermediate approximations of `U` after each iteration. return_residuals If `True`, return a |VectorArray| of all residual vectors which have been computed during the Newton iterations. Returns ------- U |VectorArray| of length 1 containing the computed solution data Dict containing the following fields: :error_sequence: |NumPy array| containing the residual norms after each iteration. :stages: See `return_stages`. :residuals: See `return_residuals`. Raises ------ NewtonError Raised if the Netwon algorithm failed to converge. """ logger = getLogger('pymor.algorithms.newton') data = {} if initial_guess is None: initial_guess = operator.source.zeros() if return_stages: data['stages'] = operator.source.empty() if return_residuals: data['residuals'] = operator.range.empty() U = initial_guess.copy() residual = rhs - operator.apply(U, mu=mu) err = residual.l2_norm()[0] if error_norm is None else error_norm(residual)[0] logger.info(f' Initial Residual: {err:5e}') iteration = 0 error_sequence = [err] while True: if iteration >= miniter: if err <= atol: logger.info(f'Absolute limit of {atol} reached. Stopping.') break if err/error_sequence[0] <= rtol: logger.info(f'Prescribed total reduction of {rtol} reached. Stopping.') break if (len(error_sequence) >= stagnation_window + 1 and err/max(error_sequence[-stagnation_window - 1:]) >= stagnation_threshold): logger.info(f'Error is stagnating (threshold: {stagnation_threshold:5e}, window: {stagnation_window}). ' f'Stopping.') break if iteration >= maxiter: raise NewtonError('Failed to converge') if iteration > 0 and return_stages: data['stages'].append(U) if return_residuals: data['residuals'].append(residual) iteration += 1 jacobian = operator.jacobian(U, mu=mu) try: correction = jacobian.apply_inverse(residual, least_squares=least_squares) except InversionError: raise NewtonError('Could not invert jacobian') U += correction residual = rhs - operator.apply(U, mu=mu) err = residual.l2_norm()[0] if error_norm is None else error_norm(residual)[0] logger.info(f'Iteration {iteration:2}: Residual: {err:5e}, ' f'Reduction: {err / error_sequence[-1]:5e}, Total Reduction: {err / error_sequence[0]:5e}') error_sequence.append(err) if not np.isfinite(err): raise NewtonError('Failed to converge') data['error_sequence'] = np.array(error_sequence) return U, data