c demonstration program for ode/de/step package. c c reference:shampine and gordon's computer solution of c ordinary differential equations:the initial value problem. c c the package is used to solve the defining equations for the c jacobian elliptic functions: c y1'=y2*y3, y1(0)=0 c y2'=-y1*y3, y2(0)=1 c y3'=-ksq*y1*y2, y3(0)=1 c which is solved for ti=i*5 for i=0, 1, ..., 12. c the analytic solution is c y1=sn(t/ksq) c y2=cn(t/ksq) c y3=dn(t/ksq) c c in this case we use ksq=k*k=0.51. c output is written on logical unit lout which has been set to 6. c real y, t, tout, relerr, abserr, work dimension y(3), work(163), iwork(5) external f data lout/6/ c neqn=3 relerr=1.0e-9 abserr=1.0e-16 iflag=1 write(lout, 1000)neqn, relerr, abserr, iflag c t=0.0e0 y(1)=0.0e0 y(2)=1.0e0 y(3)=1.0e0 write(lout, 1002)t, y, iflag do 11 i=1,12 tout=5.0e0*i 3 call ode(f, neqn, y, t, tout, relerr, abserr, iflag, work, * iwork) write(lout, 1002)t, y, iflag goto (900,11,5,7,9,900),iflag 5 write(lout, 1004) goto 3 7 write(lout, 1006) goto 900 9 write(lout, 1008) goto 900 11 continue 900 stop 1000 format('- neqn=',i3,' relerr=',1pe12.2,' abserr=',f12.2, * ' iflag=',i3,//) 1002 format(' ',4(1pe16.8,2x),i3) 1004 format('0tolerances too much and have been changed') 1006 format('0too many steps...eqn. is hard') 1008 format('0eqn. appears to be stiff') end subroutine f(t, y, yp) dimension y(3), yp(3) real t, y, yp yp(1)=y(2)*y(3) yp(2)=-y(1)*y(3) yp(3)=-0.51e0*y(1)*y(2) return end