example4_7.java

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

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

public class Example4_7 extends Applet
{
 public  static void F(int n,double[] a,double[] b,double x,double y,double[] s)
 {int i;double temp;
  s[0]=s[2]=a[0];s[1]=s[3]=b[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+b[i];
   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+b[n];
  }

 public static boolean Cssrt(int n,double[] a,double[] b,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&&b[0]==0)return false;
   //if(n==1){p=a[0]*a[0]+b[0]*b[0];x[0]=-(a[0]*a[1]+b[0]*b[1])/p;y[0]=(a[1]*b[0]-a[0]*b[1])/p;return true;}
   p=a[n]*a[n]+b[n]*b[n];
   while(p<1e-16){x[n-1]=y[n-1]=0;n--;p=a[n]*a[n]+b[n]*b[n];}
   p/=(a[0]*a[0]+b[0]*b[0]);
   u=Math.exp(Math.log(p)/2/n);p=u;
   for(i=n-1;i>=0;i--){a[i]*=p;b[i]*=p;p*=u;}
   p=a[0]*a[0]+b[0]*b[0];
   for(i=n;i>0;i--){q=(a[i]*a[0]+b[i]*b[0])/p;b[i]=(a[0]*b[i]-b[0]*a[i])/p;a[i]=q;}
   a[0]=1;b[0]=0;
   while(n>0)
   {F(n,a,b,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,b,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,b,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;}
          }
         }
    }
    x[n-1]=x0*u;y[n-1]=y0*u;
    x[n-1]=Math.round(x[n-1]*1e7)/1e7;y[n-1]=Math.round(y[n-1]*1e7)/1e7;
    for(i=1;i<n;i++){a[i]+=x0*a[i-1]-y0*b[i-1];b[i]+=x0*b[i-1]+y0*a[i-1];}
    a[n]=0;b[n]=0;
    n-=1;
    x0=y0=0;}
   return true;
 }

 public void paint(Graphics g)
 {double[] a={1,5,-90,-80,727,-645,-1262};
  double[] b={0,-15,-37,320,-136,-1265,941};
  double[] x=new double[6];
  double[] y=new double[6];
  int i;
  Cssrt(6,a,b,x,y);
  for(i=0;i<6;i++)
  {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);
  }
 }
}