example4_6.java

清华大学2002年出版的《科学与工程数值计算算法Java》配套源码

import java.applet.*;
import java.awt.*;

public class Example4_6 extends Applet
{
 public  static void F(int n,double[] a,double x,double y,double[] s)
 {int i;double temp;
  s[0]=a[0];s[1]=0;s[2]=a[0];s[3]=0;
  for(i=1;i<n;i++)
  {temp=s[0];s[0]=x*temp-y*s[1]+a[i];s[1]=x*s[1]+y*temp;
   temp=s[2];s[2]=x*temp-y*s[3]+s[0];s[3]=x*temp+y*s[3]+s[1];}
   temp=s[0];s[0]=x*temp-y*s[1]+a[n];s[1]=x*s[1]+y*temp;
  }

 public static boolean ssrt(int n,double[] a,double[] x,double[] y)
 { double[] s=new double[4];
   double x0,y0,x1,y1,p,q,p1,q1,t=1,c,d,u,dc;
   int i=0;boolean falg;
   x0=y0=x1=y1=p=q=p1=q1=0;
   if(a[0]==0)return false;
   if(n==1){x[0]=-a[1]/a[0];y[0]=0;return true;}
   if(n==2)
   {p=a[1]*a[1]-4*a[0]*a[2];
    if(p>=0){x[0]=(-a[1]+Math.sqrt(p))/2/a[0];x[1]=(-a[1]-Math.sqrt(p))/2/a[0];y[0]=y[1]=0;}
    else {x[0]=x[1]=-a[1]/2/a[0];y[0]=Math.sqrt(-p)/2/a[0];y[1]=-y[0];}
    return true;
	}
   while(Math.abs(a[n])<1e-12){x[n-1]=y[n-1]=0;n--;}
   u=Math.exp(Math.log(Math.abs(a[n]/a[0]))/n);p=u;
   for(i=n-1;i>=0;i--){a[i]*=p;p*=u;}
   for(i=n;i>0;i--)a[i]/=a[0];
   a[0]=1;
   while(n>0)
   {F(n,a,x0,y0,s);
    p=s[0]*s[0]+s[1]*s[1];
    q=s[2]*s[2]+s[3]*s[3];
    falg=true;
    while(true)
    { if(p<1e-16)break;
      if(q>1e-16&&falg)
      { t=1;
		while(true)
        {x1=x0-t*(s[0]*s[2]+s[1]*s[3])/q;
         y1=y0-t*(s[1]*s[2]-s[0]*s[3])/q;
         F(n,a,x1,y1,s);
         p1=s[0]*s[0]+s[1]*s[1];
         q1=s[2]*s[2]+s[3]*s[3];
         if(p1<p){p=p1;q=q1;x0=x1;y0=y1;t=1;break;}
         t/=2;
         if(t<1e-10){falg=false;break;}
         }
      }
     else
      { d=0.01;
		 while(true)
          { c=0;dc=Math.PI/50;
			for(i=0;i<100;i++)
           {x1=x0+d*Math.cos(c);y1=y0+d*Math.sin(c);c+=dc;
            if(x1*x1+y1*y1>1)continue;
            F(n,a,x1,y1,s);
            p1=s[0]*s[0]+s[1]*s[1];
            q1=s[2]*s[2]+s[3]*s[3];
            if(p1<p)break;}
            if(i==100)d+=0.01;
            else {p=p1;q=q1;t=1;d=0.1;x0=x1;y0=y1;falg=true;break;}
          }
         }
    }

   if(Math.abs(y0)<1e-7)
   { x[n-1]=x0*u;y[n-1]=0;
     for(i=1;i<n;i++)a[i]+=x0*a[i-1];a[n]=0;
     n-=1;
   }
   else
   { x[n-1]=x[n-2]=x0*u;y[n-1]=y0*u;y[n-2]=-y[n-1];
     p=-2*x0;q=x0*x0+y0*y0;
     for(i=1;i<n-1;i++)
     {a[i]-=a[i-1]*p;a[i+1]-=a[i-1]*q;}
     a[n-1]=a[n]=0;
     n-=2;
   }
   x0=y0=0;}
   return true;
 }

 public void paint(Graphics g)
 {double[] a={1,-13,86,-312,623,-635,250};
  double[] x=new double[6];
  double[] y=new double[6];
  int i;
  ssrt(6,a,x,y);
  for(i=0;i<6;i++)
  {x[i]=Math.round(x[i]*1e8)/1e8;y[i]=Math.round(y[i]*1e8)/1e8;
   if(y[i]==0)g.drawString("第"+(i+1)+"个根： "+x[i],10,20+20*i);
   else if(y[i]<0)g.drawString("第"+(i+1)+"个根： "+x[i]+""+y[i]+"i",10,20+20*i);
   else g.drawString("第"+(i+1)+"个根： "+x[i]+"+"+y[i]+"i",10,20+20*i);
  }
 }
}