//============================================================================ // Copyright (c) Kitware, Inc. // All rights reserved. // See LICENSE.txt for details. // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notice for more information. // // Copyright 2015 National Technology & Engineering Solutions of Sandia, LLC (NTESS). // Copyright 2015 UT-Battelle, LLC. // Copyright 2015 Los Alamos National Security. // // Under the terms of Contract DE-NA0003525 with NTESS, // the U.S. Government retains certain rights in this software. // // Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National // Laboratory (LANL), the U.S. Government retains certain rights in // this software. //============================================================================ #ifndef vtk_m_NewtonsMethod_h #define vtk_m_NewtonsMethod_h #include #include namespace vtkm { template struct NewtonsMethodResult { bool Valid; bool Converged; vtkm::Vec Solution; }; /// Uses Newton's method (a.k.a. Newton-Raphson method) to solve a nonlinear /// system of equations. This function assumes that the number of variables /// equals the number of equations. Newton's method operates on an iterative /// evaluate and search. Evaluations are performed using the functors passed /// into the NewtonsMethod. The first functor returns the NxN matrix of the /// Jacobian at a given input point. The second functor returns the N tuple /// that is the function evaluation at the given input point. The input point /// that evaluates to the desired output, or the closest point found, is /// returned. /// VTKM_SUPPRESS_EXEC_WARNINGS template VTKM_EXEC_CONT NewtonsMethodResult NewtonsMethod( JacobianFunctor jacobianEvaluator, FunctionFunctor functionEvaluator, vtkm::Vec desiredFunctionOutput, vtkm::Vec initialGuess = vtkm::Vec(ScalarType(0)), ScalarType convergeDifference = ScalarType(1e-3), vtkm::IdComponent maxIterations = 10) { using VectorType = vtkm::Vec; using MatrixType = vtkm::Matrix; VectorType x = initialGuess; bool valid = false; bool converged = false; for (vtkm::IdComponent iteration = 0; !converged && (iteration < maxIterations); iteration++) { // For Newton's method, we solve the linear system // // Jacobian x deltaX = currentFunctionOutput - desiredFunctionOutput // // The subtraction on the right side simply makes the target of the solve // at zero, which is what Newton's method solves for. The deltaX tells us // where to move to to solve for a linear system, which we assume will be // closer for our nonlinear system. MatrixType jacobian = jacobianEvaluator(x); VectorType currentFunctionOutput = functionEvaluator(x); VectorType deltaX = vtkm::SolveLinearSystem(jacobian, currentFunctionOutput - desiredFunctionOutput, valid); if (!valid) { break; } x = x - deltaX; converged = true; for (vtkm::IdComponent index = 0; index < Size; index++) { converged &= (vtkm::Abs(deltaX[index]) < convergeDifference); } } // Not checking whether converged. return { valid, converged, x }; } } // namespace vtkm #endif //vtk_m_NewtonsMethod_h