Sunday, December 9, 2012

Newton-Raphson Method, Part 2

Newton-Raphson Method, Part 2

Now we apply the procedure that we detailed in Part 1 to solve the equation:

f(x) = x^2 - 4*sin(x) = 0

First we write our own code for the functions sin() and cos(), since we are not sure of the accuracy of the built-in functions that are supplied with the math library.  We call our own versions mySin() and myCos(). Then we write the code for f() and call it fun(), and the code for f'() and call it deriv().

Here are the codes for mySin(), myCos(), fun(), and deriv().

*** funderiv.c ***
#define ABS(x) ((x) < 0.0? -(x):(x))
#define EPS  1.0e-12

// x - x^3/3! + x^5/5! - x^7/7! + ...

long double mySin(long double x)
{
  long double num = x;
  long double den = 1.0;
  long double term = x;
  long double sum = x;
  int j;
  for(j = 3; j < 1000; j += 2) {
    num = -num * x * x;
    den = den * (j-1) * j;
    term = num / den;
    sum = sum + term;
    if(ABS(term) < EPS) break;
  }
  return sum;
}

// 1 - x^2/2! + x^4/4! - x^6/6! + ...

long double myCos(long double x)
{
  long double num = 1.0;
  long double den = 1.0;
  long double term = 1.0;
  long double sum = 1.0;
  int j;
  for(j = 2; j < 1000; j += 2) {
    num = -num * x * x;
    den = den * (j-1) * j;
    term = num / den;
    sum = sum + term;
    if(ABS(term) < EPS) break;
  }
  return sum;
}

// f(x) =x^2 - 4.0*sin(x)

long double fun(long double x)
{
  return x*x - 4.0*mySin(x);
}

// f'(x) =2*x - 4*cos(x)

long double deriv(long double x)
{
  return 2.0*x - 4.0*myCos(x);
}

*** end ***

Compiling and running by using the shell script nrsolve.sh, we get the numerical results:

./nrsolve.sh funderiv.c 3.0
Starting value supplied: 3.000000
Iter=0, x=3.00000000000000, f(x)=8.43551996776053
Iter=1, x=2.15305769201338, f(x)=1.29477250528655
Iter=2, x=1.95403864200580, f(x)=0.10843855339463
Iter=3, x=1.93397153275207, f(x)=0.00115163152386
Iter=4, x=1.93375378855763, f(x)=0.00000013605495
Iter=5, x=1.93375376282702, f(x)=0.00000000000000
Solution: 1.93375376282702, error 0.00000000000100 after 5 iters


I hope that you find the procedure outlined here useful in your own work.  Enjoy!

No comments:

Post a Comment