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