My code compiles and is working correctly, but I was just wondering if anyone could check over it and suggest any improvements or correct things I've done that are bad in the c++ world...anything at all?
Thanks in advance!
#include<iostream>
#include<vector>
#include<cmath>
#include<iomanip>
using namespace std;
//Declare functions
void poly(vector<double>& A, vector<double>& B, vector<double>& C, double x, double& p, double& px, double& pxx, int& degree);
void bisection(vector<double> B, double left_endpoint, double right_endpoint, double epsilon, double p, int max, int degree, int result);
double f(vector<double>, double, double, int);
//---------------------------------------------------------------------------------------
int main()
{
vector<double> A; // vector of coefficients a^degree...a^0
vector<double> B; // backup copy of vector A
vector<double> C;
int degree; // highest degree
double x; // value of x
double p; // p(x)
double px; // p'(x)
double pxx; // p''(x)
int max = 0; // maximum number of iterations
double left_endpoint = 0.0; // left endpoint of interval
double right_endpoint = 0.0; // right endpoint of interval
double epsilon = 0.000001; // convergence criterion
int result;
poly(A, B, C, x, p, px, pxx, degree);
bisection(B, left_endpoint, right_endpoint, epsilon, p, max, degree, result);
return 0;
}
//---------------------------------------------------------------------------------------
void poly(vector<double>& A, vector<double>& B, vector<double>& C, double x, double& p, double& px, double& pxx, int& degree)
{
cout << "\nFor the polynomial p(x) = a0 + a1x + a2x^2 + ... + anx^n" << endl;
cout << "Enter the polynomial degree, n" << endl;
cout << "Example: 2 for a quadratic, 3 for a cubic..." << endl;
cout << " " << endl;
cout << "Degree: ";
cin >> degree;
A.resize (degree + 1);
B.resize (degree + 1);
C.resize (degree + 1);
for (int i = degree; i >= 0; i--)
{
cout << endl << "Enter coefficient a" << i << ": ";
cin >> A[i];
}
for (int i = degree; i >= 0; i--)
{
C[i] = B[i] = A[i];
}
cout << "\n********** DIFFERENTIATION **********"<< endl;
cout << endl <<"Enter the value of x for which to solve: ";
cin >> x;
//Calculations for p(x)
p = A[degree];
for (int i = (degree - 1); i >= 0; i--)
{
p = p*x;
p = p+A[i];
}
//First differentiation on the coefficients
for (int i = 1; i <= degree; i++)
{
A[i - 1] = i * A[i];
}
//Horner's Method - calculations for p'(x)
px = A[degree - 1];
for (int i = (degree - 2); i >= 0; i--)
{
px = px*x;
px = px+A[i];
}
//Second differentiation on the coefficients
for (int i = 1; i <= degree-1 ; i++)
{
A[i - 1] = i * A[i];
}
//Horner's Method - calculations for p''(x)
pxx = A[degree - 2];
for (int i = (degree - 3); i >= 0; i--)
{
pxx = pxx*x;
pxx = pxx+A[i];
}
cout << " p\(" << x << ") = " << p << endl;
cout << " p\'(" << x << ") = " << px << endl;
cout << "p\''(" << x << ") = " << pxx << endl;
//Reset vector A by copying values from vector B
for (int i = degree; i >= 0; i--)
{
A[i] = B[i];
}
}
//---------------------------------------------------------------------------------------
void bisection(vector<double> B, double left_endpoint, double right_endpoint, double epsilon, double p, int max, int degree, int result)
{
double x1 = 0.0; // left endpoint of current interval
double x2 = 0.0; // right endpoint of current interval
double xmid = 0.0; // computed midpoint of current interval
double f1 = 0.0; // function evaluated at left endpoint of current interval
double f2 = 0.0; // function evaluated at right endpoint of current interval
double fmid = 0.0;; // function evaluated at computed midpoint of current interval
double width = 0.0; // width of original interval (b - a)
double current_width = 0.0; // width of current interval (xmid - x1)
result = 0;
cout << "\n********** BISECTION METHOD **********"<< endl;
cout << endl << "Enter the maximum number of iterations allowed: ";
cin >> max;
for (double j = -50; j <= 50.1; j=j+0.1)
{
f(B, j, p, degree);
if ( f(B, j, p, degree) <= 0 && f(B, j+0.1, p, degree) > 0 )
{
result = 1;
cout << "\nA root is between the interval (" << j <<","<< j+0.1 <<")"<<endl;
left_endpoint = j;
x1 = left_endpoint;
f1 = f(B, x1, p, degree);
right_endpoint = j+0.1;
xmid = right_endpoint;
fmid = f(B, xmid, p, degree);
width = (right_endpoint - left_endpoint);
for (int i = 1; i <= max; i++)
{
x2 = (x1 + xmid) / 2.0;
cout << "The next midpoint, approximation is " <<fixed<<setprecision(7)<< x2 << endl;
f2 = f(B, x2, p, degree);
if (f1 * f2 <= 0.0) // root is in left half interval
{
current_width = (x2 - x1) / 2.0;
fmid = f2;
xmid = x2;
}
else // root is in right half interval
{
current_width = (xmid - x2) / 2.0;
f1 = f2;
x1 = x2;
}
if (current_width < epsilon)
{
cout << "\nA root at x = " << x2 << " was found "
<< "in " << i << " iterations" << endl;
f2 = f(B, x2, p, degree);
cout << "The value of the function at the root is " << f2 << endl;
break;
}
}
}
else if ( f(B, j, p, degree) >= 0 && f(B, j+0.1, p, degree) < 0 )
{
result = 1;
cout << "\nA root is between the interval (" << j <<","<< j+0.1 <<")"<<endl;
left_endpoint = j;
x1 = left_endpoint;
f1 = f(B, x1, p, degree);
right_endpoint = j+0.1;
xmid = right_endpoint;
fmid = f(B, xmid, p, degree);
width = (right_endpoint - left_endpoint);
for (int i = 1; i <= max; i++)
{
x2 = (x1 + xmid) / 2.0;
cout << "The next midpoint, approximation is " <<fixed<<setprecision(7)<< x2 << endl;
f2 = f(B, x2, p, degree);
if (f1 * f2 <= 0.0) // root is in left half interval
{
current_width = (x2 - x1) / 2.0;
fmid = f2;
xmid = x2;
}
else // root is in right half interval
{
current_width = (xmid - x2) / 2.0;
f1 = f2;
x1 = x2;
}
if (current_width < epsilon)
{
cout << "\nA root at x = " << x2 << " was found "
<< "in " << i << " iterations" << endl;
f2 = f(B, x2, p, degree);
cout << "The value of the function at the root is " << f2 << endl;
break;
}
}
}
}
if (result == 0)
{
cout << "\nThe search for a root has failed. The polynomial you have entered"
<< "\ncontains one or more complex roots." << endl;
}
return;
}// end of bisection
//---------------------------------------------------------------------------------------
double f(vector<double> B, double x, double p, int degree)// function to evaluate p(x)
{
p = B[degree];
for (int i = (degree - 1); i >= 0; i--)
{
p = p*x;
p = p+B[i];
}
return p;
}// end of f
//---------------------------------------------------------------------------------------