zroots.cpp

Visual C++ 常用数值算法集 源代码

void value1(double x[3],double b[3],double c[3],double dum[3],double ad[3][102])
{
	int i,j;
	for(i=1; i<=3; i++)
	{
	   x[i]=0.0;
	   b[i]=0.0;
	   c[i]=0.0;
	   dum[i]=0.0;
	}
	for(i=1; i<=3; i++)
	{
		for(j=1; j<=102; j++)
		{
			ad[i][j]=0.0;
		}
	}
}

void zroots(double a[3][6], int m, double roots[3][5], int& polish)
{
	int t,i,j,jj,m1;
	double ad[3][102],x[3],b[3],c[3],dum[3],eps,dum1;
    eps = 0.000001;
    t=0;
	m1=-1;
    for (j = 1; j<= m + 1; j++)
	{
        ad[1][j] = a[1][j];
        ad[2][j] = a[2][j];
    }
    for (j = m ;j>= 1 ;j--)
	{
        x[1] = 0;
        x[2] = 0;
        laguer(ad, j, x, eps, t);
        if (fabs(x[2]) <= 2.0 * eps *eps * fabs(x[1]))  x[2] = 0.0;
        roots[1][j] = x[1];
        roots[2][j] = x[2];
        b[1] = ad[1][j + 1];
        b[2] = ad[2][j + 1];
        for (jj = j; jj>=1; jj--)
		{
            c[1] = ad[1][jj];
            c[2] = ad[2][jj];
            ad[1][jj] = b[1];
            ad[2][jj] = b[2];
            dum1 = b[1];
            b[1] = x[1] * dum1 - x[2] * b[2] + c[1];
            b[2] = x[2] * dum1 + x[1] * b[2] + c[2];
        }
    }
    if (polish) 
	{
        for (j = 1; j<= m; j++)
		{
            dum[1] = roots[1][j];
            dum[2] = roots[2][j];
            laguer(ad, m, dum, eps, m1);
        }
    }
    for (j = 2; j<= m; j++)
	{
        x[1] = roots[1][j];
        x[2] = roots[2][j];
        for (i = j - 1; j>= 1; j--)
		{
            if (roots[1][i] <= x[1])  return;
            roots[1][i + 1] = roots[1][i];
            roots[2][i + 1] = roots[2][i];
        }
        if (roots[1][i] > x[1] ) i = 0;
        roots[1][i + 1] = x[1];
        roots[2][i + 1] = x[2];
    }
    value1(dum, c, b, x, ad);
}