Chap 9: Root finding


amirvahid
10-21-2009, 08:05 PM
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.

DavidB
10-21-2009, 09:26 PM
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.)

amirvahid
10-21-2009, 10:36 PM
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.

amirvahid
10-22-2009, 05:14 PM
Dear Dave,
What do you think?

davekw7x
10-22-2009, 05:48 PM
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

DavidB
10-22-2009, 10:05 PM
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?

amirvahid
10-22-2009, 11:55 PM
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.