C++ Newton Method

[Catch_22]

Hello,

I'm trying to solve f(x)=2-(x-1)(4-x)exp(-x)=1.8 where x is between 2.4 and 2.5.

So I started by defining g(x)=0.2-(x-1)(4-x)exp(-x) and finding where this is zero.

I decided to use Newton's method(where you can choose where you start):

/*

Calculates a solution of f(x)-1.8=0 using Newton's method.

*/

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

int main() {

double x , y;

cout << "Enter a value for x: ";
cin >> x;

for(x ; x != y ; ) {

y = x;

x = y - (0.2-(y-1)*(4-y)*exp(-y))/(-exp(-x)-3*x*exp(-x)+x*x*exp(-x));

cout << y << '\n';

}

cout << "A solution to f(x)=1.8 is " << y << '\n';

return 0;

}

A printout looks like this:

Blank@Blank-Blank ~/Documents $ ./newton
Enter a value for x: 2.45
2.45
2.47989
2.54033
2.67338
3.03523
5.49596
1.13457
1.21041
1.23695
1.24774
1.25234
1.25433
1.25521
1.25559
1.25576
1.25583
1.25587
1.25588
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
1.25589
A solution to f(x)=1.8 is 1.25589

This value works if you put it back into the equation. Out of the two possible solutions, this is the wrong one!

I tried changing the starting value but kept ending up back at this number or getting "nan" an infinite number of times.

Strangely, Octave can solve g(x)=1.8 quite easily:

octave:1> function y=f(x)
> y=0.2-(x-1).*(4-x).*exp(-x)
> endfunction
octave:2> fsolve('f', 2.5)
y = 0.015309
y = 0.015312
y = 0.015306
y = -3.9657e-05
y = -3.6945e-05
y = -4.2368e-05
y = 1.8005e-10
ans = 2.4173
octave:3>

which is the answer I want to get (here I put f(x)=g(x) as defined at the start of this post).

So, my question is: how do I get 2.4173 using Newton-Rhapson method instead of 1.25589?

Thanks in advance!

[histrungalot]

Your derivative is wrong:

// Calculates a solution of f(x)-1.8=0 using Newton's method.
#include <iostream>
#include<cmath>
using namespace std;

int main() {
   double x , y;
   cout << "Enter a value for x: ";
   cin >> x;

   for(x ; x != y ; ) {
      y = x;
      //x = y - (0.2-(y-1)*(4-y)*exp(-y))/(-exp(-x)-3*x*exp(-x)+x*x*exp(-x));
      // Your derivative is wrong
      x = y - (0.2-(x-1)*(4-x)*exp(-x))/(-(4-x)*exp(-x)+(x-1)*exp(-x)+(x-1)*(4-x)*exp(-x));
      cout << y << '\n';
   }
   cout << "A solution to f(x)=1.8 is " << y << '\n';
   return 0;
}

Output:

 ./a.out
Enter a value for x: 2.5
2.5
2.41711
2.41733
2.41733
A solution to f(x)=1.8 is 2.41733

[Catch_22]

Thanks. I thought it was a programming error. I never thought my maths would be wrong! =O