d8r1.cpp

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

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

void jacobi(double a[4][4], int n, double d[], double v[4][4], int& nrot)
{
    double b[100], z[100],g,s,c,h,tau,sss,ddd,t,tresh,sm,theta;
	int ip,iq,i,j;
    for (ip = 1;ip<=n;ip++)
	{
        for (iq = 1;iq<=n;iq++)
		{
            v[ip][iq] = 0.0;
        }
        v[ip][ip] = 1.0;
    }
    for (ip = 1; ip<=n; ip++)
	{
        b[ip] = a[ip][ip];
        d[ip] = b[ip];
        z[ip] = 0.0;
    }
    nrot = 0;
    for (i = 1;i<=50;i++)
	{
        sm = 0.0;
        for (ip = 1; ip<=n - 1; ip++)
		{
            for (iq = ip + 1; iq<=n; iq++)
			{
                sm = sm + fabs(a[ip][iq]);
            }
        }
        if (sm == 0.0)
		{
			return;
		}
        if (i < 4)
		{
            tresh = 0.2 * sm /double(n*n);
		}
        else
		{
            tresh = 0.0;
        }
        for (ip = 1; ip<=n - 1; ip++)
		{
            for (iq = ip + 1; iq<=n; iq++)
			{
                g = 100.0 * fabs(a[ip][iq]);
                sss = fabs(d[ip]) + g;
                ddd = fabs(d[iq]) + g;
                if ((i > 4) && (sss == fabs(d[ip])) && (ddd == fabs(d[iq])))
				{
                    a[ip][iq] = 0.0;
				}
                else 
				{
					if (fabs(a[ip][iq]) > tresh)
					{
						h = d[iq] - d[ip];
						if ((fabs(h) + g) == fabs(h))
						{
                        t = a[ip][iq] / h;
						}
					
                        else
						{
                        theta = 0.5 * h / a[ip][iq];
                        t = 1.0 / (fabs(theta) + sqrt(1.0 + pow(theta , 2)));
						if (theta < 0.0 ) 
						t = -t;
						}
							c = 1.0 / sqrt(1.0 + t *t);
							s = t * c;
							tau = s / (1.0 + c);
							h = t * a[ip][iq];
							z[ip] = z[ip] - h;
							z[iq] = z[iq] + h;
							d[ip] = d[ip] - h;
							d[iq] = d[iq] + h;
							a[ip][iq] = 0.0;
							for (j = 1;j<=ip - 1;j++)
								{
								g = a[j][ip];
								h = a[j][iq];
								a[j][ip] = g - s * (h + g * tau);
								a[j][iq] = h + s * (g - h * tau);
								}
							for (j = ip + 1 ;j<=iq - 1;j++)
								{
								g = a[ip][j];
								h = a[j][iq];
								a[ip][j] = g - s * (h + g * tau);
								a[j][iq] = h + s * (g - h * tau);
								}
							for (j = iq + 1 ;j<=n;j++)
								{
								g = a[ip][j];
								h = a[iq][j];
								a[ip][j] = g - s * (h + g * tau);
								a[iq][j] = h + s * (g - h * tau);
								}
							 for (j = 1 ;j<=n;j++)
								{
								g = v[j][ip];
								h = v[j][iq];
								v[j][ip] = g - s * (h + g * tau);
								v[j][iq] = h + s * (g - h * tau);
								}
							nrot = nrot + 1;
						}
				}
			}
        }
        for (ip = 1;ip<=n;ip++)
		{
            b[ip] = b[ip] + z[ip];
            d[ip] = b[ip];
            z[ip] = 0.0;
        }
    }
    cout<< "50 iterations should never happen"<<endl;
}

void main()
{
    //program d8r1
    //driver for routine jacobi
	int np,ii,jj,j,k,l,kk,ll,nrot;
    double  d[4],v[4][4],r[4];
    double  a[4][4],e[4][4],ratio;
	np = 3;
    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 (ii = 1;ii<=3;ii++)
	{
        for (jj = 1;jj<=3;jj++)
		{
            e[ii][ jj] = a[ii][ jj];
        }
    }
    jacobi(e, np,d ,v, nrot);
    cout<< "Number of jacobi rotations;  "<< nrot<<endl;
    cout<< "Eigenvalues;"<<endl;
    for (j = 1;j<=np;j++)
	{
        cout<<setw(13)<<d[j];
    }
    cout<<endl;
    cout<< "Eigenvectors;"<<endl;
    for (j = 1;j<=np;j++)
	{
        cout<< setw(5)<< "Number"<<j;
        for (k = 1;k<=np;k++)
		{
            cout<< setw(11)<<v[k][j];
        }
        cout<<endl;
    }
    cout<<endl;
    //eigenvector test
    cout<< setw(5)<< "Eigenvector test";
    for (j = 1;j<=np;j++)
	{
        for (l = 1;l<=np;l++)
		{
            r[l] = 0.0;
            for (k = 1;k<=np;k++)
			{
               if (k > l)
			   {
                    kk = l;
                    ll = k;
			   }
                else
				{
                    kk = k;
                    ll = l;
                }
                r[l] = r[l] + a[ll][kk] * v[k][j];
            }
        }
        cout<<endl;
        cout<< setw(5)<<"Vector number"<<j<<endl;
        cout<<endl;
        cout<< setw(15)<<"         vector        mtrx*vec.        ratio"<<endl;
        for (l = 1;l<=np;l++)
		{
            ratio = double(r[l]) / v[l][j];
            cout<< setw(15)<<v[l][j];
            cout<< setw(15)<<r[l];
            cout<< setw(15)<<ratio<<endl;
        }
    }
}