home *** CD-ROM | disk | FTP | other *** search
/ QBasic & Borland Pascal & C / Delphi5.iso / C / Samples / C-SSP.ARJ / DETOHA.C < prev    next >
Encoding:
Text File  |  1984-08-07  |  1.6 KB  |  51 lines
  1.    detoha(iter,h,x,y,dy,ddy,iflag,d)
  2.  
  3.       /* subroutine to solve third-order ordinary differential  */
  4.       /* equations y''' = f(x,y,y',y'') by Hamming's method.    */
  5.  
  6.       int iflag,iter;
  7.       float d[],dy[],ddy[],h,x[],y[];
  8.  
  9.    {
  10.       int i,lda;
  11.       float u2,u3,u4,pmc1,pmc2,pmc3,ppp,pmm,pp,pm,p,pmmm,c,cc,ccc,tm;
  12.       double ytrpri();
  13.  
  14.       d[0] = 0.;
  15.       d[1] = 0.;
  16.       d[2] = 0.;
  17.       d[3] = 0.;
  18.       lda = 4;
  19.       detork(lda,h,x,y,dy,ddy,iflag);
  20.       u2 = ytrpri(x[1],y[1],dy[1],ddy[1],iflag);
  21.       u3 = ytrpri(x[2],y[2],dy[2],ddy[2],iflag);
  22.       u4 = ytrpri(x[3],y[3],dy[3],ddy[3],iflag);
  23.       pmc1 = 0.;
  24.       pmc2 = 0.;
  25.       pmc3 = 0.;
  26.  
  27.       for(i = 3; i <= iter-2; i++)
  28.       {
  29.        ppp = ddy[i-3] + (4./3.) * h * (2.*u4-u3+2.*u2);
  30.        pmm =  ppp-(112./121.) * pmc1;
  31.        pp = dy[i-3] + (4./3.) * h * (2.*ddy[i]-ddy[i-1]+2.*ddy[i-2]);
  32.        pm = pp-(112./121.) * pmc2;
  33.        p = y[i-3] + (4./3.) * h * (2.*dy[i]-dy[i-1]+2.*dy[i-2]);
  34.        tm = p-(112./121.) * pmc3;
  35.        pmmm = ytrpri(x[i+1],tm,pm,pmm,iflag);
  36.        ccc = (9.*ddy[i]-ddy[i-2]+3.*h*(pmmm+2.*u4-u3))/8.;
  37.        pmc1 = ppp - ccc;
  38.        ddy[i+1] = ccc + (9./121.) * pmc1;
  39.        cc = (9.*dy[i]-dy[i-2]+3.*h*(ddy[i+1]+2.*ddy[i]-ddy[i-1]))/8.;
  40.        pmc2 = pp - cc;
  41.        dy[i+1] = cc + (9./121.) * pmc2;
  42.        c = (9.*y[i]-y[i-2]+3.*h*(dy[i+1]+2.*dy[i]-dy[i-1]))/8.;
  43.        pmc3 = p - c;
  44.        y[i+1] = c + (9./121.) * pmc3;
  45.        u2 = u3;
  46.        u3 = u4;
  47.        u4 = ytrpri(x[i+1],y[i+1],dy[i+1],ddy[i+1],iflag);
  48.        d[i+1] = c - p;
  49.       }
  50.     }
  51.