Back to index

wims  3.65+svn20090927
Quartic.java
Go to the documentation of this file.
00001 /*
00002  * Created on 18.07.2006
00003  *
00004  */
00005 package rene.zirkel.expression;
00006 
00007 public class Quartic 
00008 {
00009 
00010        static public int solve (double c[], double sol[])
00011        {      double sum=0;
00012               for (int k=0; k<c.length; k++) sum+=Math.abs(c[k]);
00013               if (sum<1e-15)
00014               {      sol[0]=0; return 1;
00015               }
00016               for (int k=0; k<c.length; k++) c[k]/=sum;
00017               int n=4;
00018               while (Math.abs(c[n])<1e-15) n--;
00019               if (n==4)
00020               {      double soli[]=new double[4];
00021                      return quartic(c,sol,soli);
00022               }
00023               else if (n==3)
00024               {      return cubic(c,sol);
00025               }
00026               else if (n==2)
00027               {      double a=-c[1]/(2*c[2]);
00028                      double r=a*a-c[0]/c[2];
00029                      if (r<-1e-10) return 0;
00030                      else if (Math.abs(r)<1e-10)
00031                      {      sol[0]=a; return 1;
00032                      }
00033                      else
00034                      {      sol[0]=a+Math.sqrt(r);
00035                             sol[1]=a-Math.sqrt(r);
00036                             return 2;
00037                      }
00038               }
00039               else if (n==1)
00040               {      sol[0]=-c[0]/c[1];
00041                      return 1;
00042               }
00043               else
00044               {      if (Math.abs(c[0])<1e-10)
00045                      {      sol[0]=0; return 1;
00046                      }
00047               }
00048               return 0;
00049        }
00050        
00051        /*-------------------- Global Function Description Block ----------------------
00052         *
00053         *     ***QUARTIC************************************************25.03.98
00054         *     Solution of a quartic equation
00055         *     ref.: J. E. Hacke, Amer. Math. Monthly, Vol. 48, 327-328, (1941)
00056         *     NO WARRANTY, ALWAYS TEST THIS SUBROUTINE AFTER DOWNLOADING
00057         *     ******************************************************************
00058         *     dd(0:4)     (i)  vector containing the polynomial coefficients
00059         *     sol(1:4)    (o)  results, real part
00060         *     soli(1:4)   (o)  results, imaginary part
00061         *     Nsol        (o)  number of real solutions 
00062         *     ==================================================================
00063         *     17-Oct-2004 / Raoul Rausch
00064         *            Conversion from Fortran to C
00065         *
00066         *
00067         *-----------------------------------------------------------------------------
00068         */
00069         static int quartic (double dd[], double sol[], double soli[])
00070         {
00071               double AA[] = new double[4], z[] = new double[3];
00072               double a, b, c, d, f, p, q, r, zsol, xK2, xL, xK, sqp, sqm;
00073               int ncube, i;
00074               int Nsol = 0;
00075 
00076               if (dd[4] == 0.0)
00077               {      return 0;
00078               }
00079 
00080               a = dd[4];
00081               b = dd[3];
00082               c = dd[2];
00083               d = dd[1];
00084               f = dd[0];
00085 
00086               p = (-3.0*Math.pow(b,2) + 8.0 *a*c)/(8.0*Math.pow(a,2));
00087               q = (Math.pow(b,3) - 4.0*a*b*c + 8.0 *d*Math.pow(a,2)) / (8.0*Math.pow(a,3));
00088               r = (-3.0*Math.pow(b,4) + 16.0 *a*Math.pow(b,2)*c - 64.0 *Math.pow(a,2)*b*d + 256.0 *Math.pow(a,3)*f)
00089                      /(256.0*Math.pow(a,4));
00090               
00091               // Solve cubic resolvent
00092               AA[3] = 8.0;
00093               AA[2] = -4.0*p;
00094               AA[1] = -8.0*r;
00095               AA[0] = 4.0*p*r - Math.pow(q,2);
00096 
00097               ncube = cubic(AA, z);
00098               
00099               zsol = - 1.e99;
00100               for (i=0;i<ncube;i++)       zsol = Math.max(zsol, z[i]);       //Not sure C has max fct
00101               z[0] =zsol;
00102               xK2 = 2.0*z[0] -p;
00103               xK = Math.sqrt(xK2);
00104               xL = q/(2.0*xK);
00105               sqp = xK2 - 4.0 * (z[0] + xL);
00106               sqm = xK2 - 4.0 * (z[0] - xL);
00107 
00108               for (i=0;i<4;i++)    soli[i] = 0.0;
00109               if ( (sqp >= 0.0) && (sqm >= 0.0))
00110               {
00111                      sol[0] = 0.5 * (xK + Math.sqrt(sqp));
00112                      sol[1] = 0.5 * (xK - Math.sqrt(sqp));
00113                      sol[2] = 0.5 * (-xK + Math.sqrt(sqm));
00114                      sol[3] = 0.5 * (-xK - Math.sqrt(sqm));
00115                      Nsol = 4;
00116               }
00117               else if ( (sqp >= 0.0) && (sqm < 0.0))
00118               {
00119                      sol[0] = 0.5 * (xK + Math.sqrt(sqp));
00120                      sol[1] = 0.5 * (xK - Math.sqrt(sqp));
00121                      sol[2] = -0.5 * xK;
00122                      sol[3] = -0.5 * xK;
00123                      soli[2] =  Math.sqrt(-.25 * sqm);
00124                      soli[3] = -Math.sqrt(-.25 * sqm);
00125                      Nsol = 2;
00126               }
00127               else if ( (sqp < 0.0) && (sqm >= 0.0))
00128               {
00129                      sol[0] = 0.5 * (-xK + Math.sqrt(sqm));
00130                      sol[1] = 0.5 * (-xK - Math.sqrt(sqm));
00131                      sol[2] = 0.5 * xK;
00132                      sol[3] = 0.5 * xK;
00133                      soli[2] =  Math.sqrt(-0.25 * sqp);
00134                      soli[3] = -Math.sqrt(-0.25 * sqp);
00135                      Nsol = 2;
00136               }
00137               else if ( (sqp < 0.0) && (sqm < 0.0))
00138               {
00139                      sol[0] = -0.5 * xK;
00140                      sol[1] = -0.5 * xK;
00141                      soli[0] =  Math.sqrt(-0.25 * sqm);
00142                      soli[1] = -Math.sqrt(-0.25 * sqm);
00143                      sol[2] = 0.5 * xK;
00144                      sol[3] = 0.5 * xK;
00145                      soli[2] =  Math.sqrt(-0.25 * sqp);
00146                      soli[3] = -Math.sqrt(-0.25 * sqp);
00147                      Nsol = 0;
00148               }
00149               
00150               for (i=0;i<4;i++)    sol[i] -= b/(4.0*a);
00151               return Nsol;
00152 
00153         }
00154         
00155         /*-------------------- Global Function Description Block ----------------------
00156          *
00157          *     ***CUBIC************************************************08.11.1986
00158          *     Solution of a cubic equation
00159          *     Equations of lesser degree are solved by the appropriate formulas.
00160          *     The solutions are arranged in ascending order.
00161          *     NO WARRANTY, ALWAYS TEST THIS SUBROUTINE AFTER DOWNLOADING
00162          *     ******************************************************************
00163          *     A(0:3)      (i)  vector containing the polynomial coefficients
00164          *     X(1:L)      (o)  results
00165          *     L           (o)  number of valid solutions (beginning with X(1))
00166          *     ==================================================================
00167          *    17-Oct-2004 / Raoul Rausch
00168          *           Conversion from Fortran to C
00169          *
00170          *
00171          *-----------------------------------------------------------------------------
00172          */
00173        static public int cubic (double A[], double X[])
00174        {
00175               final double PI = 3.1415926535897932;
00176               final double THIRD = 1./3.;
00177               double U[]=new double[3], W, P, Q, DIS, PHI;
00178               int i;
00179               int L;
00180 
00181               if (A[3] != 0.0)
00182               {
00183                      //cubic problem
00184                      W = A[2]/A[3]*THIRD;
00185                      P = Math.pow((A[1]/A[3]*THIRD - Math.pow(W,2)),3);
00186                      Q = -.5*(2.0*Math.pow(W,3)-(A[1]*W-A[0])/A[3] );
00187                      DIS = Math.pow(Q,2)+P;
00188                      if ( DIS < 0.0 )
00189                      {
00190                             //three real solutions!
00191                             //Confine the argument of ACOS to the interval [-1;1]!
00192                             PHI = Math.acos(Math.min(1.0,Math.max(-1.0,Q/Math.sqrt(-P))));
00193                             P=2.0*Math.pow((-P),(5.e-1*THIRD));
00194                             for (i=0;i<3;i++)    U[i] = P*Math.cos((PHI+2*((double)i)*PI)*THIRD)-W;
00195                             X[0] = Math.min(U[0], Math.min(U[1], U[2]));
00196                             X[1] = Math.max(Math.min(U[0], U[1]),Math.max( Math.min(U[0], U[2]), Math.min(U[1], U[2])));
00197                             X[2] = Math.max(U[0], Math.max(U[1], U[2]));
00198                             L = 3;
00199                      }
00200                      else
00201                      {
00202                             // only one real solution!
00203                             DIS = Math.sqrt(DIS);
00204                             X[0] = CBRT(Q+DIS)+CBRT(Q-DIS)-W;
00205                             L=1;
00206                      }
00207               }
00208               else if (A[2] != 0.0)
00209               {
00210                      // quadratic problem
00211                      P = 0.5*A[1]/A[2];
00212                      DIS = Math.pow(P,2)-A[0]/A[2];
00213                      if (DIS > 0.0)
00214                      {
00215                             // 2 real solutions
00216                             X[0] = -P - Math.sqrt(DIS);
00217                             X[1] = -P + Math.sqrt(DIS);
00218                             L=2;
00219                      }
00220                      else
00221                      {
00222                             // no real solution
00223                             L=0;
00224                      }
00225               }
00226               else if (A[1] != 0.0)
00227               {
00228                      //linear equation
00229                      X[0] =A[0]/A[1];
00230                      L=1;
00231               }
00232               else
00233               {
00234                      //no equation
00235                      L=0;
00236               }
00237         /*
00238          *     ==== perform one step of a newton iteration in order to minimize
00239          *          round-off errors ====
00240          */
00241               for (i=0;i<L;i++)
00242               {
00243                      X[i] = X[i] - (A[0]+X[i]*(A[1]+X[i]*(A[2]+X[i]*A[3])))/(A[1]+X[i]*(2.0*A[2]+X[i]*3.0*A[3]));
00244               }
00245 
00246               return L;
00247        }
00248 
00249        static double CBRT (double Z)
00250        {
00251               final double THIRD = 1./3.;
00252               if (Z>0) return Math.pow(Z,THIRD);
00253               else if (Z<0) return -Math.pow(Math.abs(Z),THIRD);
00254               else return 0;
00255        }
00256 
00260        public static void main(String[] args) 
00261        {      double a[]={0,3,-3,-6,5},x[]=new double[5],y[]=new double[5];
00262               int n=quartic(a,x,y);
00263               System.out.println(n+" solutions!");
00264               for (int i=0; i<n; i++) System.out.println(x[i]);
00265        }
00266 
00267 }