c This program solves the equations for phi and delta phi using an analytic c approximation for - actually two such approximations one with a c positive and one with a negative sign. c IMPLICIT NONE Integer it,nt,JUNK, nbad,nok Real*8 t,ti,Y,phiplus,phiminus,phi0,delphiplus,delphiminus,delphiexplus Real*8 delphiexminus,errphiclass Real*8 delphi0,delphicplus,delphicminus,delphicexplus,delphicexminus Real*8 eps,h1,hmin Real*8 delphiex,delphicex,delt Real*8 phiclass,phi1plus,phi1minus,difplusdelphi,difplusdelphic Real*8 difminusdelphi,difminusdelphic DIMENSION Y(14) External DERIVS, BSSTEP OPEN(UNIT = 11, FILE = 'delphian.dat') Open(unit=12,file='delphiplus.out') Open(unit=13,file='delphiminus.out') READ(11,*) JUNK READ(11,*) delt,nt,phi0,delphi0 READ(11,*) JUNK Read(11,*) eps,h1,hmin c We will start at t = 0 c Note that the equation for delphian is linear so it really doesn't matter c what the initial value is. c c I will set it up so that the first thing I do is integrate the harmonic c oscillator as a check: so that is Y(1) and Y(2). Then I'll do the c equation for phi with 1+phi^2 and then delphi. Then I'll do the equation c for phi with 1-phi^2 and then delphi. Y(1)= 1.d0 Y(2)= 0.d0 Y(3) = phi0 Y(4) = 0.d0 Y(5)= delphi0 Y(6)= 0.d0 Y(7) = phi0 Y(8) = 0.d0 Y(9)= delphi0 Y(10)= 0.d0 c Now do the exact ones for a larger phi0 Y(11) = phi0+delphi0 Y(12) = 0.d0 Y(13) = phi0+delphi0 Y(14) = 0.d0 phiclass = Y(1) phiplus = Y(3) delphiplus = Y(5) phiminus = Y(7) delphiplus = Y(9) errphiclass = dcos(t)-Y(1) phi1plus = Y(11) phi1minus = Y(13) delphicplus = delphiplus-(delphi0/phi0)*phiplus delphicminus = delphiminus-(delphi0/phi0)*phiminus delphiexplus = phi1plus - phiplus delphiexminus = phi1minus - phiminus delphicexplus = delphiexplus-(delphi0/phi0)*phiplus difplusdelphi = delphiexplus-delphiplus difplusdelphic = delphicexplus-delphicplus difminusdelphi = delphiexminus-delphiminus difminusdelphic = delphicexminus-delphicminus t = 0.d0 WRITE(12,35) t,errphiclass,phiplus,difplusdelphi,delphiplus, 1 difplusdelphic,delphicplus WRITE(13,35) t,errphiclass,phiminus,difminusdelphi,delphiminus, 1 difminusdelphic,delphicminus c I set t = 0 above DO IT = 1,nt ti = t t = t+delt CALL ODEINT(Y,14,ti,t,EPS,H1,HMIN,NOK,NBAD,DERIVS,BSSTEP) phiclass = Y(1) phiplus = Y(3) delphiplus = Y(5) phiminus = Y(7) delphiminus = Y(9) errphiclass = dcos(t)-Y(1) phi1plus = Y(11) phi1minus = Y(13) delphicplus = delphiplus-(delphi0/phi0)*phiplus delphicminus = delphiminus-(delphi0/phi0)*phiminus delphiexplus = phi1plus - phiplus delphiexminus = phi1minus - phiminus delphicexplus = delphiexplus-(delphi0/phi0)*phiplus delphicexminus = delphiexminus-(delphi0/phi0)*phiminus difplusdelphi = delphiexplus-delphiplus difplusdelphic = delphicexplus-delphicplus difminusdelphi = delphiexminus-delphiminus difminusdelphic = delphicexminus-delphicminus WRITE(12,35) t,errphiclass,phiplus,difplusdelphi,delphiplus, 1 difplusdelphic,delphicplus WRITE(13,35) t,errphiclass,phiminus,difminusdelphi,delphiminus, 1 difminusdelphic,delphicminus Enddo 35 Format(1x,10(d14.7,1x)) STOP END ! -------------------------------------------------------------- ! -------------- SUBROUTINES ----------------------------------- ! -------------------------------------------------------------- SUBROUTINE DERIVS(t,Y,DYDT,NVAR) Implicit None INTEGER nvar REAL*8 t, Y(nvar), DYDT(nvar) Real*8 phiplus,delphiplus,phiminus,delphiminus,phi1plus,phi1minus phiplus = Y(3) delphiplus = Y(5) phiminus = Y(7) delphiminus = Y(9) phi1plus = Y(11) phi1minus = Y(13) DYDT(1) = Y(2) DYDT(2) = -Y(1) DYDT(3) = Y(4) DYDT(4) = -(1.d0+phiplus**2)*phiplus DYDT(5) = Y(6) DYDT(6) = -(1.d0+3.d0*phiplus**2)*delphiplus DYDT(7) = Y(8) DYDT(8) = -(1.d0-phiminus**2)*phiminus DYDT(9) = Y(10) DYDT(10) = -(1.d0-3.d0*phiminus**2)*delphiminus DYDT(11) = Y(12) DYDT(12) = -(1.d0+phi1plus**2)*phi1plus DYDT(13) = Y(14) DYDT(14) = -(1.d0-phi1minus**2)*phi1minus RETURN END