Thursday, June 24, 2010

Using Newton method find root - C++/Matlab

1) C++
--------------------------------------------------------------------------------
/* An example of Newton's method
Author: Jiansen Lu
Created: June, 24, 2010
___________________________________

Purpose: This program estimates the root of the function f(x)


Algorithm:   The improvements to our initial guess x0  are obtained by
repeated  use of the formula:

x_new = x_old - ( f(x_old) / f'(x_old) ).

The initial guess ( x0 ), the maximum number of iterations(nmax),
the root tolerance (xtol), the function value tolerance (ftol)
are all inputed through the keyboard.

The root tolerance is defined as follows:
The iterations will stop if

| x_new - x_old | / | x_old | < xtol.


The function value tolerance is defined as follows:
The iterations will stop if

| f( x_new) | < ftol.

For accuracy all variables will be declared as double.
*/

#include <iostream>
#include <iomanip>
#include <cmath>
 using namespace std;

double fun(double x);
double funprime(double x);

int main()
{
double x_old = 0.6, x_new = 0 ;  //  estimates of the root

 double xtol = 1e-5 , ftol = 1e-5 ;  // tolerances

 double f = 0. , fprime = 0. ;   // values of the function and of its derivative

 int nmax = 1000 ;                  // max number of iterations
 int n = 0 ;                     // running index



// Prepare the heading of the table


cout <<"\n\n    x_old          f(x_old)       f'(x_old)           x_new   \n\n";
cout<< setiosflags(ios:: scientific);

cout<< setprecision(10);

//    START THE LOOP

for( n = 1; n <= nmax ; n++ )
{

f = fun(x_old);
fprime = funprime(x_old);

x_new = x_old - f / fprime ;

cout<<setw(20)<< x_old << setw(20) << f << setw(20) << fprime;

cout<<setw(20)<< x_new << endl;

if( fabs((x_new - x_old ) / x_old ) < xtol )   break;

if( fabs(f) < ftol )  break;

x_old = x_new ;
}
cout<<endl<<endl;
return (0);
}

double fun(double x){
return 2*x*x -x;
}

double funprime(double x){
return 4*x-1;
}
-----------------------------------------------------------------
2) Matlab:
f = @(x)2*x.^2-x;
Then find the zero near 0.6:
z = fzero(f,0.6)

No comments:

Post a Comment