Hi,
I am using the code zrhqr.h in order to obtain all of the roots for a 4-th order polynomial. After some MD simulations, it gives me the following error:
ERROR: Too many iterations in hqr
in file Ewald.cpp at line 201
terminate called after throwing an instance of 'int'
Aborted

Here is my code:

void PredictEvent (Int na, Int nb)//cf. Rapaport chap14
{
VecR dr, dv, rs, shift, tm, w, da;
VecI m1v, m2v;
Doub b, d, t, tInt, vv, aa, ar, av, minrr = 0;
Int dir, evCode, m1x, m1y, m1z, n;
const int pol = 4;
const int poll = pol + 1;
VecComplex rt(pol);
VecDoub Coeffs(poll);
Int indxmin = -1;
VCopy (w, mol[na].inCell);
if (mol[na].rv.x > 0.) ++ w.x;
if (mol[na].rv.y > 0.) ++ w.y;
if (mol[na].rv.z > 0.) ++ w.z;
VMul (w, w, region);
VDiv (w, w, cells);
VSAdd (rs, mol[na].r, 0.5, region);
VVSub (w, rs);
WhenCross (x);
WhenCross (y);
WhenCross (z);
if (tm.y < tm.z) dir = (tm.x < tm.y) ? 0 : 1;
else dir = (tm.x < tm.z) ? 0 : 2;
evCode = 100 + dir;
ScheduleEvent (na, MOL_LIMIT + evCode, timeNow + VComp (tm, dir));
for (m1z = cellRange[0].z; m1z <= cellRange[1].z; m1z ++) {
m1v.z = m1z;
VWrapEvC (z);
for (m1y = cellRange[0].y; m1y <= cellRange[1].y; m1y ++) {
m1v.y = m1y;
VWrapEvC (y);
for (m1x = cellRange[0].x; m1x <= cellRange[1].x; m1x ++) {
m1v.x = m1x;
VWrapEvC (x);
n = VLinear (m2v, cells) + nMol;
for (n = cellList[n]; n >= 0; n = cellList[n]) {
if (n != na && n != nb && (nb >= -1 || n < na)) {
tInt = timeNow - mol[n].time;
VSub (dr, mol[na].r, mol[n].r);
VVSAdd (dr, - tInt, mol[n].rv);
VVSAdd (dr, -0.5 * Sqr (tInt), mol[n].ra);
VVSub (dr, shift);
VSub (dv, mol[na].rv, mol[n].rv);
VVSAdd (dv, - tInt, mol[n].ra);
VSub (da, mol[na].ra, mol[n].ra);
b = VDot (dr, dv);
vv = VLenSq (dv);
aa = VLenSq (da);
ar = VDot (da, dr);
av = VDot (da, dv);
Coeffs[0] =VLenSq(dr)-1;
Coeffs[1] = 2*b;
Coeffs[2] = vv + ar;
Coeffs[3] = av;
Coeffs[4] = 0.25 * aa;

if (b < 0. && abs(Coeffs[4]) < 1e-8) {
vv = VLenSq (dv);
d = Sqr (b) - vv * (VLenSq (dr) - 1.);
if (d >= 0.) {
t = - (sqrt (d) + b) / vv;
ScheduleEvent (na, n, timeNow + t);
}
}

if (b < 0. && abs(Coeffs[4])>1e-8) {
zrhqr(Coeffs,rt);
for (int m1 = 0; m1 < pol; m1++){
if (imag(rt[m1])== 0 && real(rt[m1]) > 0) {
minrr = real (rt[m1]);
indxmin = m1;
//break;
}
if (imag(rt[m1]) == 0 && real(rt[m1]) < minrr && real(rt[m1]) > 0){
minrr = real(rt[m1]);
indxmin = m1;
}
}
if (indxmin != -1){

t = minrr;
ScheduleEvent (na, n, timeNow + t);
}
}
}
}
}
}
}
}



I marked it with "red".

Here is a general code that one can use for finding polynomial root using numerical recipes
#include "nr3.h"
#include "eigen_unsym.h"
#include "zrhqr.h"
int main()
{
const int n = 6; // degree of polynomial
Doub coeff[] = { 10.0, -18.0, 25.0, -24.0, 16.0, -6.0, 1.0};
const int N = n+1;
VecDoub a(N, coeff);
//cout << a[2] << endl;
VecComplex rt(n);
cout << endl << "Roots of polynomial x^6-6x^5+16x^4-24x^3+25x^2-18x+10" << endl;
zrhqr(a,rt);
cout << fixed << setprecision(6);

for (int i = 0; i < n; i++){
cout << setw(5) << noshowpos << i;
cout << setw(25) << showpos << rt[i] << endl;
}

return 0;
}

I'll get the following results:

Roots of polynomial x^6-6x^5+16x^4-24x^3+25x^2-18x+10
0 (+2.000000,+1.000000)
1 (+2.000000,-1.000000)
2 (+1.000000,+1.000000)
3 (+1.000000,-1.000000)
4 (+0.000000,+1.000000)
5 (+0.000000,-1.000000)


I wonder if you could help me.

I am curious: did any of your simulations run successfully?

Until you get your code figured out, perhaps an online solver would get you out of a bind:

Polynomial Root-finder (http://www.akiti.ca/PolyRootRe.html)

(At least you can get some numbers, and not have your work halted, while you check your code.)

Yes. All of my simulations run correctly when the # of charges are small in the system. I think I should somehow modify the numerical recipes code. I wonder if you could help me based on the error message that I got:

ERROR: Too many iterations in hqr
in file Ewald.cpp at line 201
terminate called after throwing an instance of 'int'
Aborted

Thanks.

Dear Dave,
What do you think?

Dear Dave,
What do you think?

I'm not an expert on zrhqr.

Here's a thought: Print out the values of the coefficients of the problematic polynomial. Try the rootfinder suggest by DavidB. What does it show you?

Ill-conditioned polynomials crop up in the most inconvenient places and I think many libraries have several different routines to try to find roots if other methods fail. (The Jenkins-Traub method, for example, is mentioned in the NR book.)

Regards,

Dave

Yes. All of my simulations run correctly when the # of charges are small in the system.

Your first post stated that you are seeking the roots for a 4-th order polynomial. However, your sample problem is a 6-th degree polynomial, for which you seem to have found the correct roots. Is this a different problem?

What is the data you enter that produces the error message?

Oh. That was just a typical example that I have extracted from NR example book. However, my simulations are using 4-th order polynomials. I increased the number of iterations "its" in "eigen_unsym.h". Right now everything seems alright. We'll see what happens next. Thanks guys for your support and elaborations.