home *** CD-ROM | disk | FTP | other *** search
/ QBasic & Borland Pascal & C / Delphi5.iso / C / Samples / C-SSP.ARJ / NAREGU.C < prev    next >
Encoding:
Text File  |  1984-07-31  |  1.7 KB  |  91 lines
  1.    naregu(xlow,xhi,eps,root,n)
  2.  
  3.       /*this subroutine finds the root of the equation f(x) = 0  */
  4.       /*through the method of false position(regula falsi).      */
  5.  
  6.       int *n;
  7.       float eps,*root,xhi,xlow;
  8.  
  9.    {
  10.       int nmax;
  11.       float fx,fxlow,fxhi,fxnew,fxp,fxq,x,xnew,xp,xq;
  12.       extern double fabs();
  13.  
  14.       nmax = *n;
  15.       *n = 0;
  16.       value(xhi,&fxhi);
  17.       value(xlow,&fxlow);
  18.       if(fxhi < 0.0)
  19.       {
  20.       xp = xlow;
  21.       fxp = fxlow;
  22.       xq = xhi;
  23.       fxq = fxhi;
  24.       }
  25.       else
  26.        {
  27.        xp = xhi;
  28.        fxp = fxhi;
  29.        xq = xlow;
  30.        fxq = fxlow;
  31.        }
  32. reit: xnew = xp - fxp*(xp - xq)/(fxp - fxq);
  33.       value(xnew,&fxnew);
  34.       *n = *n + 1;
  35.       if(fabs(fxnew) <= eps)
  36.       {
  37.       *root = xnew;
  38.       *n = 1;
  39.       return;
  40.       }
  41.       if(*n >= nmax)
  42.       {
  43.       *root = xnew;
  44.       *n = 2;
  45.       return;
  46.       }
  47.       if(fxnew >  0.0)
  48.       {
  49.       x = xnew - fxnew*(xnew - xp)/(fxnew - fxp);
  50.       xp = xnew;
  51.       fxp = fxnew;
  52.       }
  53.        else
  54.        {
  55.         x = xnew - fxnew*(xnew-xq)/(fxnew - fxq);
  56.         xq = xnew;
  57.         fxq = fxnew;
  58.        }
  59.       if(xp <= xq)
  60.       {
  61.        if(xq <= x || x <= xp) goto reit;
  62.       }
  63.       else if(xp <= x || x <= xq) goto reit;
  64.  
  65.       *n = *n + 1;
  66.       value(x,&fx);
  67.       if(fabs(fx) <= eps)
  68.       {
  69.       *root = x;
  70.       *n = 1;
  71.       return;
  72.       }
  73.       if(*n >= nmax)
  74.       {
  75.       *root = x;
  76.       *n = 2;
  77.       return;
  78.       }
  79.       if(fx <= 0.0)
  80.       {
  81.       xq = x;
  82.       fxq = fx;
  83.       }
  84.        else
  85.        {
  86.         xp = x;
  87.         fxp = fx;
  88.        }
  89.       goto reit;
  90.    }
  91.