yakabi.txt

雅克比算法

#include <stdio.h>
#include <math.h>
#include <iostream.h>
#include <fstream.h>
//#include <iomanip.h>

#define n 3
#define z 0.0001


int e,f;


int main(void)
{ 
   
  int i,j,k;  
  double a[n][n]={{3,-4,3},{-4,6,3},{3,3,1}};
  double y[n][n]={{0,0,0},{0,0,0},{0,0,0}};
  double max,p=z+1;
  double d,t,c,s;
  double maximum(double b[n][n]);
  double rest(double c[n][n]);
  while(p>z)
  {  
     max=maximum(a);// find out the max
     //cout<<max<<endl;
     //if(max=0)
       //cout<<"error"; 
     d=(a[e][e]-a[f][f])/(2*max);
     if(d<0)
       t=(-1/(fabs(d)+sqrt(d*d+1)));
     else
       t=1/(d+sqrt(d*d+1));
     c=1/(sqrt(1+t*t));
     s=c*t;
     y[e][e]=a[e][e]+t*a[e][f];
     y[f][f]=a[f][f]-t*a[e][f];
     y[e][f]=y[f][e]=0;
     //a[e][e]=a[e][e]*c*c+a[f][f]*s*s+2*max*c*s;
     //a[f][f]=a[e][e]*s*s+a[f][f]*c*c-2*max*c*s;
     //a[e][f]=a[f][e]=0;
     for(i=0;i<n;i++)
       if((i!=e)&(i!=f))
       {
           y[e][i]=a[e][i]*c+a[f][i]*s;
           y[i][e]=a[i][e]*c+a[i][f]*s;
       }            
     for(j=0;j<n;j++)
       if((j!=e)&(j!=f))
       {
           y[f][j]=a[f][j]*c-a[e][j]*s; 
           y[j][f]=a[j][f]*c-a[j][e]*s;
           
       }
     for(i=0;i<n;i++)
         for(j=0;j<n;j++)
           if((i!=e)&(i!=f)&(j!=e)&(j!=f))
               y[i][j]=a[i][j];
     for(i=0;i<n;i++)
         for(j=0;j<n;j++)
           a[i][j]=y[i][j];       
     p=rest(a); 
      //for(i=0;i<n;i++)
        //for(j=0;j<n;j++)
       // cout<<"a["<<i<<"]["<<j<<"]="<<a[i][j]<<endl;
  }  
  for(i=0;i<n;i++)
    cout<<"a["<<i<<"]["<<i<<"]="<<a[i][i]<<endl;  
    
    
}

double maximum(double b[n][n])
{
  int i,j;
  double q=0;
  for(i=0;i<n;i++)
     for(j=0;j<n;j++)
       if(i!=j)
         {
           if(fabs(b[i][j])>fabs(q))
             { 
                 q=b[i][j];                 
                 e=i;
                 f=j;
                 
             }    
         }
  //cout<<e<<endl;
  //cout<<f<<endl;       
  return(q);           
}    
double rest(double c[n][n])
{
  int i,j;
  double x=0;
  for(i=0;i<n;i++)
     for(j=0;j<n;j++)
       if(i!=j)
       x+=c[i][j]*c[i][j];  
  return(x); 
    
}