home *** CD-ROM | disk | FTP | other *** search
/ QBasic & Borland Pascal & C / Delphi5.iso / C / Samples / C-SSP.ARJ / DETORK.C < prev    next >
Encoding:
Text File  |  1984-08-07  |  1.4 KB  |  45 lines
  1.    detork(iter,h,x,y,dy,ddy,iflag)
  2.  
  3.       /* subroutine to solve third-order ordinary differential  */
  4.       /* equations y''' = f(x,y,y',y'') by Runge-Kutta method.  */
  5.  
  6.       int iflag,iter;
  7.       float dy[],ddy[],h,x[],y[];
  8.  
  9.     {
  10.       int i;
  11.       float dyn,ddyn,tk0,tk1,tk2,tk3,tm0,tm1,tm2,tm3,xn,yn,
  12.             tn0,tn1,tn2,tn3;
  13.       double ytrpri();
  14.  
  15.       for(i = 0; i <= iter-2; i++)
  16.       {
  17.        tk0 = h * dy[i];
  18.        tn0 = h * ddy[i];
  19.        tm0 = h * ytrpri(x[i],y[i],dy[i],ddy[i],iflag);
  20.        xn = x[i] + (h/2.);
  21.        yn = y[i] + (tk0/2.);
  22.        dyn = dy[i] + (tn0/2.);
  23.        ddyn = ddy[i] + (tm0/2.);
  24.        tk1 = h * (dy[i] + (tn0/2.));
  25.        tn1 = h * (ddy[i] + (tm0/2.));
  26.        tm1 = h * ytrpri(xn,yn,dyn,ddyn,iflag);
  27.        yn = y[i] + (tk1/2.);
  28.        dyn = dy[i] + (tn1/2.);
  29.        ddyn = ddy[i] + (tm1/2.);
  30.        tk2 = h * (dy[i] + (tn1/2.));
  31.        tn2 = h * (ddy[i] + (tm1/2.));
  32.        tm2 = h * ytrpri(xn,yn,dyn,ddyn,iflag);
  33.        xn = x[i] + h;
  34.        yn = y[i] + tk2;
  35.        dyn = dy[i] + tn2;
  36.        ddyn = ddy[i] + tm2;
  37.        tk3 = h * (dy[i] + tn2);
  38.        tn3 = h * (ddy[i] + tm2);
  39.        tm3 = h * ytrpri(xn,yn,dyn,ddyn,iflag);
  40.        ddy[i+1] = ddy[i]+(tm0+2.*tm1+2.*tm2+tm3)/6;
  41.        dy[i+1] = dy[i] + (tn0+2.*tn1+2.*tn2+tn3)/6;
  42.        y[i+1] = y[i] +(tk0+2.*tk1+2.*tk2+tk3)/6.;
  43.       }
  44.     }
  45.