home *** CD-ROM | disk | FTP | other *** search
- detoma(iter,h,x,y,dy,ddy,iflag,d)
-
- /* subroutine to solve third-order ordinary differential */
- /* equations y''' = f(x,y,y',y'') by modified Adams' method.*/
-
- int iflag,iter;
- float d[],dy[],ddy[],h,x[],y[];
-
- {
- int i,lda;
- float ddyi,dyi,u1,u2,u3,u4,yi;
- double ytrpri();
-
- d[0] = 0.;
- d[1] = 0.;
- d[2] = 0.;
- lda = 4;
- detork(lda,h,x,y,dy,ddy,iflag);
- u1 = ytrpri(x[0],y[0],dy[0],ddy[0],iflag);
- u2 = ytrpri(x[1],y[1],dy[1],ddy[1],iflag);
- u3 = ytrpri(x[2],y[2],dy[2],ddy[2],iflag);
- u4 = ytrpri(x[3],y[3],dy[3],ddy[3],iflag);
-
- for(i = 3; i <= iter-2; i++)
- {
- ddyi = ddy[i]+h*(55.*u4-59.*u3+37.*u2-9.*u1)/24.;
- dyi = dy[i]+h*(55.*ddy[i]-59.*ddy[i-1]+37.*ddy[i-2]-9.*ddy[i-3])/24.;
- yi = y[i]+h*(55.*dy[i]-59.*dy[i-1]+37.*dy[i-2]-9.*dy[i-3])/24.;
- u1 = u2;
- u2 = u3;
- u3 = u4;
- u4 = ytrpri(x[i+1],yi,dyi,ddyi,iflag);
- ddy[i+1] = ddy[i]+h*(9.*u4+19.*u3-5.*u2+u1)/24.;
- dy[i+1] = dy[i]+h*(9.*ddy[i+1]+19.*ddy[i]-5.*ddy[i-1]+ddy[i-2])/24.;
- y[i+1] =y[i]+h*(9.*dy[i+1]+19.*dy[i]-5.*dy[i-1]+dy[i-2])/24.;
- u4 = ytrpri(x[i+1],y[i+1],dy[i+1],ddy[i+1],iflag);
- d[i+1] = -19.*(y[i+1]-yi)/270.;
- }
- }