NonlinearSystem Class Reference

#include <NonlinearSystem.h>

Collaboration diagram for NonlinearSystem:

Public Member Functions
int	apply ()
vector_type	compute_at (const vector_type &x)
int	compute_jacobians (const vector_type &x, bool compute_inverse=true)
void	correction_step ()
	Makes the first guess. After this call the Newton method. More...
bool	find_roots_gnewton (double tolerance=1e-7, size_t max_iter=50)
	Find roots using Newton-Raphson method. More...
template<typename T >
void	initialize (const T &x)
	NonlinearSystem (size_t systemSize)
value_type	slope (unsigned int which_dimen)
	Compute the slope of function in given dimension. More...
int	system (const vector_type &x, vector_type &f)
string	to_string ()

Public Attributes
double	dx_
vector_type	f_
matrix_type	invJ_
bool	is_f_positive_
bool	is_jacobian_valid_
matrix_type	J_
ReacInfo	ri
const size_t	size_
vector_type	slopes_
vector_type	x1
vector_type	x2
vector_type	x_

Detailed Description

Definition at line 88 of file NonlinearSystem.h.

Constructor & Destructor Documentation

NonlinearSystem::NonlinearSystem ( size_t  systemSize )

inline

Definition at line 93 of file NonlinearSystem.h.

                                         : size_( systemSize )
    {
        f_.resize( size_, 0);
        slopes_.resize( size_, 0);
        x_.resize( size_, 0 );

        J_.resize( size_, size_, 0);
        invJ_.resize( size_, size_, 0);

        x2.resize( size_, 0);
        x1.resize( size_, 0);

        ri.nVec.resize( size_ );
        dx_ = sqrt( numeric_limits<double>::epsilon() );
    }

Member Function Documentation

int NonlinearSystem::apply ( )

inline

Definition at line 116 of file NonlinearSystem.h.

    {
        return system(x_, f_);
    }

vector_type NonlinearSystem::compute_at ( const vector_type &  x )

inline

Definition at line 109 of file NonlinearSystem.h.

    {
        vector_type result( size_ );
        system(x, result);
        return result;
    }

int NonlinearSystem::compute_jacobians ( const vector_type &  x, bool  compute_inverse = true  )

inline

Definition at line 121 of file NonlinearSystem.h.

References inverse().

    {
        for( size_t i = 0; i < size_; i++)
            for( size_t j = 0; j < size_; j++)
            {
                vector_type temp = x;
                temp[j] += dx_;
                J_(i, j) = (compute_at(temp)[i] - compute_at(x)[i]) / dx_;
            }

        // is_jacobian_valid_ = true;
        // Keep the inverted J_ ready
        //if(is_jacobian_valid_ and compute_inverse )
        if( compute_inverse )
            inverse( J_, invJ_ );

        return 0;
    }

Here is the call graph for this function:

void NonlinearSystem::correction_step ( )

inline

Makes the first guess. After this call the Newton method.

Definition at line 279 of file NonlinearSystem.h.

    {
        // Get the jacobian at current point. Notice that in this method, we
        // don't have to compute inverse of jacobian

        vector_type direction( size_ );

        // Now take the largest step possible such that the value of system at
        // (x_ - step ) is lower than the value of system as x_.
        vector_type nextState( size_ );

        apply();

        unsigned int i = 0;

        double factor = 1e-2;
        while( true )
        {
            i += 1;
            compute_jacobians( x_, false );
            // Make a move in either side of direction. In whichever direction
            // the function decreases.
            direction = ublas::prod( J_, f_ );
            nextState = x_ - factor * direction;
            if( ublas::norm_2( compute_at( nextState ) ) >= ublas::norm_2(compute_at(x_)))
                factor = factor / 2.0;
            else
            {
                cerr << "Correction term applied ";
                x_ = nextState;
                apply();
                break;
            }

            if ( i > 20 )
                break;
        }
    }

bool NonlinearSystem::find_roots_gnewton ( double  tolerance = 1e-7, size_t  max_iter = 50  )

inline

Find roots using Newton-Raphson method.

Parameters

tolerance	1e-7
max_iter	Maximum number of iteration allowed , default 100

Returns

If successful, return true. Check the variable x_ at which the system f_ is close to zero (within the tolerance).

Definition at line 225 of file NonlinearSystem.h.

    {
        //tolerance = sqrt( numeric_limits<double>::epsilon() );
        double norm2OfDiff = 1.0;
        size_t iter = 0;
        int status = apply();

        while( ublas::norm_2(f_) >= tolerance )
        {
            iter += 1;
            compute_jacobians( x_, true );
            vector_type correction = ublas::prod( invJ_, f_ );
            x_ -=  correction;

            // If could not compute the value of system successfully.
            status = apply();
            if( 0 != status )
                return false;

            if( iter >= max_iter )
                break;

        }

        ri.nIter = iter;

        if( iter >= max_iter )
            return false;

        return true;
    }

template<typename T >

void NonlinearSystem::initialize ( const T &  x )

inline

Definition at line 141 of file NonlinearSystem.h.

References init().

    {
        vector_type init;
        init.resize(size_, 0);

        for( size_t i = 0; i < size_; i++)
            init[i] = x[i];

        x_ = init;
        apply();

        compute_jacobians( init );
    }

Here is the call graph for this function:

value_type NonlinearSystem::slope ( unsigned int  which_dimen )

inline

Compute the slope of function in given dimension.

Parameters

which_dimen

The index of dimension.

Returns

Slope.

Definition at line 264 of file NonlinearSystem.h.

    {
        vector_type x = x_;
        x[which_dimen] += dx_;
        // x1 and x2 holds the f_ of system at x_ and x (which is x +
        // some step)
        system( x_, x1 );
        system( x, x2 );
        return ublas::norm_2( (x2 - x1)/dx_ );
    }

int NonlinearSystem::system ( const vector_type &  x, vector_type &  f  )

inline

Definition at line 168 of file NonlinearSystem.h.

    {
        int num_consv = ri.num_mols - ri.rank;
        for ( size_t i = 0; i < ri.num_mols; ++i )
        {
            double temp = x[i] * x[i] ;

#if 0
            // if overflow
            if ( std::isnan( temp ) or std::isinf( temp ) )
            {
                cerr << "Failed: ";
                for( auto v : ri.nVec ) cerr << v << ", ";
                cerr << endl;
                return -1;
            }
#endif
            ri.nVec[i] = temp;
        }

        vector< double > vels;
        ri.pool->updateReacVelocities( &ri.nVec[0], vels );

        assert( vels.size() == static_cast< unsigned int >( ri.num_reacs ) );

        // y = Nr . v
        // Note that Nr is row-echelon: diagonal and above.
        for ( int i = 0; i < ri.rank; ++i )
        {
            double temp = 0;
            for ( int j = i; j < ri.num_reacs; ++j )
                temp += ri.Nr(i, j ) * vels[j];
            f[i] = temp ;
        }

        // dT = gamma.S - T
        for ( int i = 0; i < num_consv; ++i )
        {
            double dT = - ri.T[i];
            for ( size_t  j = 0; j < ri.num_mols; ++j )
                dT += ri.gamma(i, j) * x[j] * x[j];

            f[ i + ri.rank] = dT ;
        }
        return 0;
    }

string NonlinearSystem::to_string ( )

inline

Definition at line 155 of file NonlinearSystem.h.

    {
        stringstream ss;

        ss << "=======================================================";
        ss << endl << setw(25) << "State of system: " ;
        ss << " Argument: " << x_ << " Value : " << f_;
        ss << endl << setw(25) << "Jacobian: " << J_;
        ss << endl << setw(25) << "Inverse Jacobian: " << invJ_;
        ss << endl;
        return ss.str();
    }

Member Data Documentation

double NonlinearSystem::dx_

Definition at line 321 of file NonlinearSystem.h.

vector_type NonlinearSystem::f_

Definition at line 323 of file NonlinearSystem.h.

matrix_type NonlinearSystem::invJ_

Definition at line 327 of file NonlinearSystem.h.

bool NonlinearSystem::is_f_positive_

Definition at line 330 of file NonlinearSystem.h.

bool NonlinearSystem::is_jacobian_valid_

Definition at line 329 of file NonlinearSystem.h.

matrix_type NonlinearSystem::J_

Definition at line 326 of file NonlinearSystem.h.

ReacInfo NonlinearSystem::ri

Definition at line 335 of file NonlinearSystem.h.

const size_t NonlinearSystem::size_

Definition at line 319 of file NonlinearSystem.h.

vector_type NonlinearSystem::slopes_

Definition at line 325 of file NonlinearSystem.h.

vector_type NonlinearSystem::x1

Definition at line 333 of file NonlinearSystem.h.

vector_type NonlinearSystem::x2

Definition at line 333 of file NonlinearSystem.h.

vector_type NonlinearSystem::x_

Definition at line 324 of file NonlinearSystem.h.

---

The documentation for this class was generated from the following file:

• moose-core/ksolve/NonlinearSystem.h