d8r4.cpp

这个是常用的数值算法用VC编写的。相信对大家有用哦。

# include <math.h>
# include <iomanip.h>
# include <iostream.h>

void tqli(double d[4],double  e[4],int  n, double z[4][4])
{
	int i,l,m,iter,k,t;
	double dd,g,r,s,c,p,f,b;
    if (n > 1)
	{
        for (i = 2; i<=n; i++)
		{
            e[i - 1] = e[i];
        }
        e[n] = 0.0;
        for (l = 1; l<=n; l++)
		{
            iter = 0;
loop1:          for (m = l; m<=n - 1; m++)
				{
                dd = fabs(d[m]) + fabs(d[m + 1]);
                if ((fabs(e[m]) + dd) == dd ) goto  loop2;
				}
				m = n;
loop2:			if (m != l) 
				{
                if (iter == 30 ) cout<< " too many iterations ";
                iter = iter + 1;
                g = (d[l + 1] - d[l]) / (2.0 * e[l]);
                r = sqrt(g * g + 1.0);
				if (g>0) t=1;
				if (g==0) t=0;
				if (g<0) t=-1;
                g = d[m] - d[l] + e[l] / (g + fabs(r) * t);
                s = 1.0;
                c = 1.0;
                p = 0.0;
                for (i = m - 1 ; i>=l; i--)
				{
                    f = s * e[i];
                    b = c * e[i];
                    if (fabs(f) >= fabs(g))
					{
                        c = g / f;
                        r = sqrt(c *c + 1.0);
                        e[i + 1] = f * r;
                        s = 1.0 / r;
                        c = c * s;
					}
                    else
					{
                        s = f / g;
                        r =sqrt(s *s + 1.0);
                        e[i + 1] = g * r;
                        c = 1.0 / r;
						s=s*c;
                    }
                    g = d[i + 1] - p;
                    r = (d[i] - g) * s + 2.0 * c * b;
                    p = s * r;
                    d[i + 1] = g + p;
                    g = c * r - b;
                    //omit lines from here ...
                    for (k = 1 ;k<= n;k++)
					{
                        f = z[k][i + 1];
                        z[k][i + 1] = s * z[k][i] + c * f;
                        z[k][i] = c * z[k][i] - s * f;
                    }
                   //to here when finding only eigenvalues.
                }
                d[l] = d[l] - p;
                e[l] = g;
                e[m] = 0.0;
                goto loop1;
				}
		}
    }
}

void tred2(double a[4][4], int n, double d[4], double e[4])
{
	int i,j,k,l;
	double h,scale1,f,g,hh,t;
    if (n > 1)
	{
        for (i = n; i>=2; i--)
		{
            l = i - 1;
            h = 0.0;
            scale1 = 0.0;
            if (l > 1)
			{
                for (k = 1; k<=l; k++)
					scale1 = scale1 + fabs(a[i][k]);
                if (scale1 == 0.0)
					e[i] = a[i][l];
			    else
				{
                    for (k = 1; k<=l; k++)
					{
                        a[i][ k] = a[i][k] / scale1;
                        h = h + a[i][k]*a[i][k];
                    }
                    f = a[i][ l];
					if (f>0.0) t=1;
					if (f==0.0) t=0;
					if (f<0.0)  t=-1;
                    g = -sqrt(h) * t;
                    e[i] = scale1 * g;
                    h = h - f * g;
                    a[i][ l] = f - g;
                    f = 0.0;
                    for (j = 1; j<=l; j++)
					{
                        a[j][i] = a[i][j] / h;
                        g = 0.0;
                        for (k = 1; k<=j; k++)
					    g = g + a[j][k] * a[i][k];
                        if (l > j )
						{
                           for (k = j + 1; k<=l; k++)
						   g = g + a[k][j] * a[i][k];
                        }
                        e[j] = g / h;
                        f = f + e[j] * a[i][j];
                    }
                    hh = f / (h + h);
                    for (j = 1 ; j<=l; j++)
					{
                        f = a[i][j];
                        g = e[j] - hh * f;
                        e[j] = g;
                        for (k = 1; k<=j; k++)
						a[j][k] = a[j][k] - f * e[k] - g * a[i][k];
                    }
               }
			}
            else
			{
                e[i] = a[i][l];
            }
            d[i] = h;
		}
    }
    //omit following line if finding only eigenvalues.
    d[1] = 0.0;
    e[1] = 0.0;
    for (i = 1; i<=n; i++)
	{
    //delete lines from here ...
        l = i - 1;
        if (d[i] != 0.0)
		{
            for (j = 1 ; j<=l; j++)
			{
                g = 0.0;
                for (k = 1 ; k<=l; k++)
				{
                    g = g + a[i][ k] * a[k][j];
                }
                for (k = 1 ;k<=l;k++)
				{
                    a[k][j] = a[k][j] - g * a[k][i];
                }
            }
		}
    //... to here when finding only eibenvalues.
        d[i] = a[i][i];
    //also delete lines from here ...
        a[i][i] = 1.0;
        if (l >= 1)
		{
            for (j = 1 ; j<=l; j++)
			{
                a[i][j] = 0.0;
                a[j][i] = 0.0;
            }
        }
    //... to here when finding only eigenvalues.
    }
}

void main()
{
    //program d8r4
    //driver for routine tqli
    int np,k,i,j;
    double tiny,a[4][4], c[4][4], d[4], e[4], f[4];;
	np = 3;
	tiny = 0.000001;
    a[1][1] = 1.0; a[1][2] = 2.0; a[1][3] = 3.0;
    a[2][1] = 2.0; a[2][2] = 2.0; a[2][3] = 3.0;
    a[3][1] = 3.0; a[3][2] = 3.0; a[3][3] = 3.0;
    for (i = 1; i<=np; i++)
	{
        for (j = 1; j<=np; j++)
		{
            c[i][j] = a[i][j];
        }
    }
    tred2(c, np, d, e);
    tqli(d, e, np, c);
    cout<<endl;
    cout<< "Eigenvectors for a real symmetric matrix";
    for (i = 1; i<=np; i++)
	{
        for (j = 1; j<=np; j++)
		{
            f[j] = 0.0;
            for (k = 1; k<=np; k++)
			{
                f[j] = f[j] + a[j][k] * c[k][i];
            }
		}
        cout<<endl;
        cout<< "Eigenvalue "<< i<< " = ";
        cout<< setw(10)<<d[i];
        cout<<endl;
        cout<<endl;
        cout<<  "       vector     mtrx*vect.    ratio"<<endl;
        cout<<endl;
        for (j = 1; j<=np; j++)
		{
            if (fabs(c[j][i]) < tiny )
			{
              cout<<setw(10)<<c[j][i];
              cout<< setw(10)<<f[j];
              cout<< "div. by 0";
			  cout<<endl;
			}
			
            else
			{
              cout<<setw(13)<<c[j][ i];
              cout<<setw(13)<<f[j];
              cout<<setw(13)<<(f[j] / c[j][ i])<<endl;
            }
		}
	
    }
}