program ode c c n=number of nonlinear odes c nn=n*(n+1)=total number of odes c nn1=nn+1 (for compilers that start arrays at element 0) c parameter (n=3, nn=12, nn1=13) c implicit double precision(a-h,o-z) c external fcn c dimension y(nn1),yprime(nn1),v(nn1),A(nn1),B(nn1),C(nn1) dimension D(nn1),cum(n),znorm(n),gsc(n) c c NSTEP is the total number of reorthonormalization steps that c will be performed. Set this to any very large value. c IRATE is the number of numerical integration time steps per c reorthonormalization. c STPSZE is the integration stepsize in seconds. Choose a # c of "seconds" small compared to the "mean orbital period". c IO is the rate at which output is generated. Rather than c outputing the spectrum for each of the NSTEP reorthonorms, c values are generated for every IOth reorthonorm. c write(*,111) 111 format(1x,'nstep, stpsze, irate, io : ') read(*,*) nstep,stpsze,irate,io c c initial conditions for nonlinear ODEs c (** Choose within the system's basin of attraction **) c v(1)=1.0 v(2)=1.0 v(3)=1.0 tme=0.0 c c initial conditions for linearized ODEs c (** Leave these #s alone! They are problem independent! **) c do 10 i=n+1,nn v(i)=0.0 10 continue do 20 i=1,n v((n+1)*i)=1.0 cum(i)=0.0 20 continue c do 100 m=1,nstep c do 25 j=1,irate c c ***************************************************** c do 26 i=1,nn y(i)=v(i) 26 continue t=tme call fcn(t, y, yprime) do 27 i=1,nn A(i)=yprime(i) 27 continue c c ***************************************************** c do 28 i=1,nn y(i)=v(i)+(stpsze*A(i))/2.0 28 continue t=tme+stpsze/2.0 call fcn(t, y, yprime) do 29 i=1,nn B(i)=yprime(i) 29 continue c c ***************************************************** c do 30 i=1,nn y(i)=v(i)+(stpsze*B(i))/2.0 30 continue t=tme+stpsze/2.0 call fcn(t, y, yprime) do 31 i=1,nn C(i)=yprime(i) 31 continue c c ***************************************************** c do 32 i=1,nn y(i)=v(i)+(stpsze*C(i)) 32 continue t=tme+stpsze call fcn(t, y, yprime) do 33 i=1,nn D(i)=yprime(i) 33 continue c c ***************************************************** c do 34 i=1,nn v(i)=v(i)+stpsze*(A(i)+D(i)+2.0*(B(i)+C(i)))/6.0 34 continue tme=tme+stpsze c c ***************************************************** c 25 continue c c construct new orthonormal basis by gram-schmidt c c normalize first vector c znorm(1)=0.0 do 38 j=1,n znorm(1)=znorm(1)+v(n*j+1)**2 38 continue znorm(1)=sqrt(znorm(1)) do 40 j=1,n v(n*j+1)=v(n*j+1)/znorm(1) 40 continue c c generate new orthonormal set c c make j-1 gsr coefficients c do 80 j=2,n c do 50 k=1,j-1 gsc(k)=0.0 do 50 l=1,n gsc(k)=gsc(k)+v(n*l+j)*v(n*l+k) 50 continue c c construct a new vector c do 60 k=1,n do 60 l=1,j-1 v(n*k+j)=v(n*k+j)-gsc(l)*v(n*k+l) 60 continue c c calculate the vector's norm c znorm(j)=0.0 do 70 k=1,n znorm(j)=znorm(j)+v(n*k+j)**2 70 continue znorm(j)=sqrt(znorm(j)) c c normalize the new vector c do 80 k=1,n v(n*k+j)=v(n*k+j)/znorm(j) 80 continue c c update running vector magnitudes c do 90 k=1,n cum(k)=cum(k)+alog(znorm(k))/alog(2.0) 90 continue c c normalize exponent and print every io iterations c if(mod(m,io).ne.0) goto 100 write(5,334) tme,(cum(k)/tme,k=1,n) 334 format(1x,f12.6,2x,e12.6,2x,e12.6,2x,e12.6) c 100 continue c end c c *************************************************************** c subroutine fcn(t,y,yprime) c implicit double precision(a-h,o-z) c dimension y(13),yprime(13) c c parameters for Lorenz ODEs c b=4.0 sg=16.0 r=45.92 c c Nonlinear Lorenz equations c yprime(1)=sg*(y(2)-y(1)) yprime(2)=-y(1)*y(3)+r*y(1)-y(2) yprime(3)=y(1)*y(2)-b*y(3) c c Linearized Lorenz equations c do 10 i=0,2 yprime(4+i)=sg*(y(7+i)-y(4+i)) yprime(7+i)=(r-y(3))*y(4+i)-y(7+i)-y(1)*y(10+i) yprime(10+i)=y(2)*y(4+i)+y(1)*y(7+i)-b*y(10+i) 10 continue c return end