Hi again,

I make minor changes to the xqroot program of NR 2.11 C++, as shown below. The quadratics of x^4 - 1 is x^2 -1 and x^2 + 1.

when I run the following program I get:
<b,c> = <0,-1> and <0,0>

You should get:
<b,c> = <0,-1> and <0,1>

Sincerely,
King


#include "stdafx.h" // default C++ compiling, VC++
#include <iostream>
#include <iomanip>
#include <complex>
#include <cmath>
#include "nr.h"
using namespace std;


int _tmain(int argc, _TCHAR* argv[])
{
const int N=4; // degree of polynomial
const int NTRY=10;
const DP EPS=1.0e-6,TINY=1.0e-5;
const DP p_d[N+1]={-1,0.0,0.0,0.0,1.0};
bool newroot;
int i,j,nroot=0;
Vec_DP p(p_d,N+1),b(NTRY),c(NTRY);

cout << endl << "P(x)=x^4 - 1" << endl;
cout << "Quadratic factors x^2+bx+c" << endl << endl;
cout << setw(6) << "factor" << setw(11) << "b";
cout << setw(13) << "c" << endl << endl;
cout << fixed << setprecision(6);
for (i=0;i<NTRY;i++) {
c[i]=0.5*(i+1);
b[i] = -0.5*(i+1);
NR::qroot(p,b[i],c[i],EPS);
if (nroot == 0) {
cout << setw(4) << nroot << setw(16) << b[i];
cout << setw(13) << c[i] << endl;
nroot=1;
} else {
newroot=true;
for (j=0;j<nroot;j++)
if ((fabs(b[i]-b[j]) < TINY) && (fabs(c[i]-c[j]) < TINY))
newroot=false;
if (newroot) {
cout << setw(4) << nroot << setw(16) << b[i];
cout << setw(13) << c[i] << endl;
++nroot;
}
}
}
return 0;
}

Dear King,

As the book explains, qroot only converges properly if the initial guess is good. You can't simply use the b[i] and c[i] from a different problem. Try b[i] = EPS. Also, you may need c[i] positive for the one root, negative for the other.

The text points out that qroot is best for polishing approximate roots found by another method.

Saul Teukolsky

Hi again,

I make minor changes to the xqroot program of NR 2.11 C++, as shown below. The quadratics of x^4 - 1 is x^2 -1 and x^2 + 1.

when I run the following program I get:
<b,c> = <0,-1> and <0,0>

You should get:
<b,c> = <0,-1> and <0,1>

Sincerely,
King


#include "stdafx.h" // default C++ compiling, VC++
#include <iostream>
#include <iomanip>
#include <complex>
#include <cmath>
#include "nr.h"
using namespace std;


int _tmain(int argc, _TCHAR* argv[])
{
const int N=4; // degree of polynomial
const int NTRY=10;
const DP EPS=1.0e-6,TINY=1.0e-5;
const DP p_d[N+1]={-1,0.0,0.0,0.0,1.0};
bool newroot;
int i,j,nroot=0;
Vec_DP p(p_d,N+1),b(NTRY),c(NTRY);

cout << endl << "P(x)=x^4 - 1" << endl;
cout << "Quadratic factors x^2+bx+c" << endl << endl;
cout << setw(6) << "factor" << setw(11) << "b";
cout << setw(13) << "c" << endl << endl;
cout << fixed << setprecision(6);
for (i=0;i<NTRY;i++) {
c[i]=0.5*(i+1);
b[i] = -0.5*(i+1);
NR::qroot(p,b[i],c[i],EPS);
if (nroot == 0) {
cout << setw(4) << nroot << setw(16) << b[i];
cout << setw(13) << c[i] << endl;
nroot=1;
} else {
newroot=true;
for (j=0;j<nroot;j++)
if ((fabs(b[i]-b[j]) < TINY) && (fabs(c[i]-c[j]) < TINY))
newroot=false;
if (newroot) {
cout << setw(4) << nroot << setw(16) << b[i];
cout << setw(13) << c[i] << endl;
++nroot;
}
}
}
return 0;
}

Hello King

A ready made solution for functional equation for finding the root etc is brought forward that might be of use under your contextual situation mentioned by you
This may not be the exact requirement for yours though epsilon etc is also included(as mentioned by the author) so please have a try
;)




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

void maximum(double, double, double, double); // function prototype
double f(double); // function prototype

int main()
{
double delta; // step size
double a, b; // left and right ends of the original interval
double epsilon; // convergence criterion

// obtain the input data
cout << "Enter the limits of the original search interval, a and b: ";
cin >> a >> b;
cout << "Enter the convergence criteria: ";
cin >> epsilon;
cout << "Enter the step size: ";
cin >> delta;

maximum(a, b, epsilon, delta);

cin.get();cin.ignore(); // this line is optional

return 0;
}

void maximum(double a, double b, double epsilon, double delta)
{
double x1, x2;
double f1, f2;
double max;

// echo back the passed input data
cout << "\nThe original search interval is from " << a << " to " << b << endl;
cout << "The convergence criterion is: interval < " << epsilon << endl;
cout << "The step size is " << delta << endl;

x1 = a;
max = x1;
while(delta >= epsilon && max < b)
{
x2 = x1 + delta;
f1 = f(x1);
f2 = f(x2);
if (f1 < f2)
{
x1 = x2;
max = x2;
}
else
{
delta = delta / 2;
max = x1;
}
}

if(max > b && max != b)
max = max - delta;
cout << "\nThe max was found at x = " << max << endl;

return;
}

// function to evaluate f(x)
double f(double x)
{
const double PI = 2*asin(1.0); // value of pi

return (exp(-x) - sin(0.5 * PI * x));
}