Shyam Guthikonda

//--------------------------------------------------------------------
// Name: Shyam M Guthikonda
// Email: shyamguth@gmail.com
// URL: http://shy.am
// Desc: Program to calculate the root of a non-linear function, utilizing
//		 the Newton-Raphson method.
//--------------------------------------------------------------------
 
//--------------------------------------------------------------------
// Include Files
//--------------------------------------------------------------------
#include "math.h"
#include "windows.h"
 
//--------------------------------------------------------------------
// Globals
//--------------------------------------------------------------------
unsigned int	Iter;			// Keeps track of the number of iterations each
								// root location method requires. Avoids the need
								// to pass an extra parameter by reference.
 
ULONG			ErrorFlags = 0;	// Used to keep track of any errors.
 
// Possible error flags
enum ERRORS
{
	ERRORS_DERIV_EQ_ZERO	= 1,
	ERRORS_DENOM_EQ_ZERO	= 2,
	ERRORS_MAX_ITERATION	= 4,
	ERROR_FORCE_32BIT		= 0x7FFFFFFF
};
 
// Struct to hold the X and Y roots of a non-linear system
struct Point
{
	double X;
	double Y;
};
 
//--------------------------------------------------------------------
// Name: F1()
// Desc: The first function.
//--------------------------------------------------------------------
double F1( double x, double y )
{
	return ( (x * x) + 1 - y );
}
 
//--------------------------------------------------------------------
// Name: F2()
// Desc: The second function.
//--------------------------------------------------------------------
double F2( double x, double y )
{
	return ( 3 * cos( x ) - y );
}
 
//--------------------------------------------------------------------
// Name: F1dx()
// Desc: Derivative of the first function with respect to x.
//--------------------------------------------------------------------
double F1dx( double x, double y )
{
	return ( 2 * x );
}
 
//--------------------------------------------------------------------
// Name: F1dy()
// Desc: Derivative of the first function with respect to y.
//--------------------------------------------------------------------
double F1dy( double x, double y )
{
	return ( -1 );
}
 
//--------------------------------------------------------------------
// Name: F2dx()
// Desc: Derivative of the second function with respect to x.
//--------------------------------------------------------------------
double F2dx( double x, double y )
{
	return ( -3 * sin( x ) );
}
 
//--------------------------------------------------------------------
// Name: F2dy()
// Desc: Derivative of the second function with respect to y.
//--------------------------------------------------------------------
double F2dy( double x, double y )
{
	return ( -1 );
}
 
//--------------------------------------------------------------------
// Name: Es()
// Desc: Calculates the percent tolerance, given a desired number of
//		 significant figures.
//--------------------------------------------------------------------
double Es( unsigned int SigFigs )
{
	return (.5 * pow( 10.0, (2.0 - (double)SigFigs) ));
}
 
//--------------------------------------------------------------------
// Name: NewtonRaphson()
// Desc: Determines the root of a non-linear function using the
//		 Newton-Raphson method.
//--------------------------------------------------------------------
Point NewtonRaphsonNonLinear( double Xi, double Yi, double Es, unsigned int Imax )
{
	double Ea = 1000.0, Denom;
	double Xrold = Xi, Yrold = Yi;
	Point Root;
	Iter = 0;
 
	Root.X = Xi;
	Root.Y = Yi;
 
	// No errors yet with this method!
	ErrorFlags = 0;
 
	while ( 1 )
	{
		Xrold = Root.X; Yrold = Root.Y;
		Iter++;
 
		// Prevent division by 0.
		Denom = F1dx( Root.X, Root.Y ) * F2dy( Root.X, Root.Y ) - F1dy( Root.X, Root.Y ) * F2dx( Root.X, Root.Y );
		if (Denom == 0.0)
		{
			// Set the given error flag and return the most recent root estimates.
			ErrorFlags |= ERRORS_DENOM_EQ_ZERO;
			return Root;
		}
 
		// Newton-Raphson formula for a non-linear system.
		Root.X = Xrold - ( (F1( Xrold, Yrold ) * F2dy( Xrold, Yrold )) - (F2( Xrold, Yrold ) * F1dy( Xrold, Yrold )) ) / Denom;
		Root.Y = Yrold - ( (F2( Xrold, Yrold ) * F1dx( Xrold, Yrold )) - (F1( Xrold, Yrold ) * F2dx( Xrold, Yrold )) ) / Denom;
 
		// Calculate the relative absolute error, Ea (of X term)
		// Don't calculate the error on the first iteration since there is no previous term
		if ( (Root.X != 0.0) && (Iter != 1) ) Ea = abs( (Root.X - Xrold) / Root.X ) * 100;
 
		// Root found within error range desired, or number of iterations has
		// reached the maximum amount permitted by Imax. Either way: return.
		if ( Ea < Es	  )	{ return Root;									   }
		if ( Iter >= Imax ) { ErrorFlags |= ERRORS_MAX_ITERATION; return Root; }
 
	} // End while
}
 
//--------------------------------------------------------------------
// Name: ErrorMsg()
// Desc: Function to print the currently set bits in ErrorFlags, of enum
//		 type ERRORS.
//--------------------------------------------------------------------
void ErrorMsg()
{
	if ( ErrorFlags )
	{
		printf( "Errors:\n" );
		if ( ErrorFlags & ERRORS_DERIV_EQ_ZERO  ) printf( "\tDerivative = 0\n"			);
		if ( ErrorFlags & ERRORS_DENOM_EQ_ZERO  ) printf( "\tDenominator = 0\n"			);
		if ( ErrorFlags & ERRORS_MAX_ITERATION  ) printf( "\tMax Iterations Reached\n"	);
	} // End if errors
}
 
//--------------------------------------------------------------------
// Name: _tmain() (Application Entry Point)
// Desc: Entry point for the program.
//--------------------------------------------------------------------
int _tmain()
{
	Point		Root;
	double		fEs = Es( 6 );	// Number of desired significant figures.
 
	printf( "\nRoots Accurate to 6 Sig Figs\n\n" );
 
	// NEWTON-RAPHSON-----------------------------------------------
	Root = NewtonRaphsonNonLinear(
					1.5,			// Initial X guess
					3.5,			// Initial Y guess
					fEs,			// Percent tolerance
					1000			// Maximum iterations allowed
					);
 
	printf( "Method: Newton-Raphson (Non-Linear)\n" );
 
	// Display any errors.
	ErrorMsg();
 
	printf( "X_Root = %g\n", Root.X );
	printf( "Y_Root = %g\n", Root.Y );
	printf( "Iterations = %i\n", Iter );
 
	return 0;
}