Newton Raphson method cubic polynomial

Sep 10, 2011 at 9:59am
Hey, I am a novice to this, and now have to build a program to read sequentially, line by line, the coefficients of a cubic polynomial, and the initial guess, approximate the root and display the number of iterations it took for it to be correct to within 5e-6, or if it fails to converge, to add 10 to the estimate and try again, if it the fails, display a fail message and move to the next line of the input file, do it all again and again until it reaches the end of input file and quits.

In short, I am more than a little screwed, and will barter my firstborn for help.

here is my code so far.

#include "stdafx.h"
#include <math.h>
#include <complex>

#define a 1 /* The coefficient a, is always one, so is a defined constant */
#define EPSILON 5.0e-6

double b, c, d, n, x, xnew, guess, poly, primepoly; /*this line declares variables */
int count, *addcount;
double *addb, *addc, *addd, *addn, *addx, *addxnew, *addguess;

int
main (void)
{
addb = &b;
addc = &c;
addd = &d;
addcount = &count;
addn = &n;
addx = &x;
addxnew = &xnew;
addcount = &count;


printf("This program calculates the root of a cubic polynomial\n");

FILE *inp;
inp = fopen("E:\\input1.txt", "r+");
scanf("%lf %lf %lf %lf"); *addb, *addc, *addd, *addguess;
{
count = *addcount ++;
poly = (a * x * x * x) + (b * x * x)+(c * x) + d; /* F(x) */
primepoly = (3 * x * x) + (b * x) + c; /* F'(x) */

*addxnew = x - poly/primepoly;
}
printf("The solution is x = %.5", xnew);
return (0);
}
Sep 10, 2011 at 1:06pm
Hello~~ A cubic polynomial has at least one real root while your solution may sometimes fail to find one.
I wrote a toy program which solves this problem, but not completely.
Here's the code. This program runs correctly in Windows.
#include <iostream>
#include <assert.h>
using namespace std;

double const EPS = 5.0e-10;
double const STEP = 10;

double a(0), b(0), c(0), d(0);
double x(0);
double root(0);
double step(10);

inline double F(double x) { return (a*x*x*x + b*x*x + c*x + d); }
inline double f(double x) { return (3*a*x*x + 2*b*x + c); }

inline bool Judge(double Fx) { return ( Fx < EPS && Fx > -EPS ); } // Fx == 0 ?

int main()
{
unsigned int count(0);
cout << "input a, b, c, d : ";
cin >> a; cin >> b; cin >> c; cin >> d;
assert(cin.good());
assert( a < -EPS || a > EPS );

cout << "initial guess : ";
cin >> x;
assert(cin.good());

double preX(0), curX(0), preF(0), curF(0);
curX = preX = x;
do
{
curF = F(curX);
if ( Judge(curF) )
{
root = curX;
break;
}

double flag = f(curX);

if (curF > 0 && flag > 0)
{
preF = curF;
step = -STEP;
preX = curX;
curX = curX + step;
curF = F(curX);
count++;
}
else if (curF > 0 && flag < 0)
{
preF = curF;
preX = curX;
curX = curX + step;
curF = F(curX);
count++;
}
else if (curF < 0 && flag > 0)
{
preF = curF;
preX = curX;
curX = curX + step;
curF = F(curX);
count++;
}
else
{
preF = curF;
step = -STEP;
preX = curX;
curX = curX + step;
curF = F(curX);
count++;
}
cout << "count = " << count << endl;
} while (preF*curF > 0);

double bigFx = curF;
double bigX = curX;
double smallFx = preF;
double smallX = preX;
if (preF > 0)
{
bigFx = preF;
bigX = preX;
smallFx = curF;
smallX = curX;
}

double mediumX = (bigX + smallX) / 2;
while ( !Judge(F(mediumX)) )
{
if ( F(mediumX) > 0 )
{
bigX = mediumX;
}
else
{
smallX = mediumX;
}
mediumX = (bigX + smallX) / 2;
double tmp = F(mediumX);
}

root = mediumX;
cout << "root = " << root << endl;
return 0;
}

Good Luck and Have a Nice Day:)
Last edited on Sep 10, 2011 at 1:06pm
Sep 10, 2011 at 3:01pm
Thanks heaps, that was really helpful to see the way you strung together the functions, I am actually getting correct answers! now i just need to get the data from the .dat file and assign the values, then loop the calculator until the end of the file.

It is just a bit much for a first year assignment for me...
Sep 10, 2011 at 3:08pm
sorry to have made a mistake in my code.
the code segment if...else if....else if...else should be replaced by the following code.
if (curF > 0)
{
preF = curF;
step = -STEP;
preX = curX;
curX = curX + step;
curF = F(curX);
count++;
}
else
{
preF = curF;
preX = curX;
curX = curX + STEP;
curF = F(curX);
count++;
}
Calculating the differentiation of F(x) is of no help.
Topic archived. No new replies allowed.