NewtonRaphson.h

Go to the documentation of this file.

#pragma once

#include <Eigen/Sparse>
#include "nuto/base/Exception.h"
#include "nuto/math/LineSearch.h"

namespace NuTo
{
namespace NewtonRaphson
{

class NoConvergence : public NuTo::Exception
{
public:
    NoConvergence(const std::string& caller = "", const std::string& message = "")
        : Exception(caller, message)
    {
    }
};

struct DoubleSolver
{
    static double Solve(double dr, double r)
    {
        return r / dr;
    }
};

template <typename TR, typename TDR, typename TNorm, typename TTol, typename TInfo>
struct Problem
{
    TR Residual;
    TDR Derivative;
    TNorm Norm;
    TTol mTolerance;
    TInfo Info;
};

template <typename TR, typename TDR, typename TNorm, typename TTol, typename TInfo = VoidInfo>
auto DefineProblem(TR residual, TDR derivative, TNorm norm, TTol tolerance, TInfo info = VoidInfo())
{
    return Problem<TR, TDR, TNorm, TTol, TInfo>({residual, derivative, norm, tolerance, info});
}

template <typename TNonlinearProblem, typename TX, typename TSolver, typename TLineSearchAlgorithm = NoLineSearch>
auto Solve(TNonlinearProblem&& problem, TX&& x0, TSolver&& solver, int maxIterations = 20,
           TLineSearchAlgorithm&& lineSearch = NoLineSearch(), int* numIterations = nullptr)
{
    auto x = x0;
    auto r = problem.Residual(x);

    int iteration = 0;
    problem.Info(iteration, x, r);

    if (problem.Norm(r) < problem.mTolerance)
        return x;

    while (iteration < maxIterations)
    {
        auto dr = problem.Derivative(x);
        auto dx = solver.Solve(dr, r);

        ++iteration;

        if (lineSearch(problem, &r, &x, dx))
        {
            if (numIterations)
                *numIterations = iteration;
            problem.Info(iteration, x, r);
            return x;
        }
        problem.Info(iteration, x, r);
    }
    if (numIterations)
        *numIterations = iteration;
    throw NoConvergence(__PRETTY_FUNCTION__, "No convergence after " + std::to_string(iteration) + " iterations.");
}
} /* NewtonRaphson */
} /* NuTo */