* * K-dV eq. with periodic boundary condition * * Using Runge-Kutta-Gill method / Ver 3.1 * * Programed by : T. Tamaribuchi * IMPLICIT REAL*8 (A-H,O-Z) REAL xmin, ymin, xmax, ymax, ch INTEGER NOG PARAMETER(KBLUE = 16, LSOLID = 1, NOTXORPEN = 11) COMMON A,B,T,ET,DT,XS,DX,ND,NPD,J COMMON /PLT/xmin, ymin, xmax, ymax, ch, NOG * CALL initial CALL GWSETBK(IA, -2) CALL GWANCHOR(IR, IA+1) * NOG = NOG + 1 CALL GWsetogn(IR, 0, -NOG) DO WHILE(INT((ET-T)*DT/DT**2) .GT. 0) CALL UPDATE IF(MOD(J,NPD) .EQ.0) THEN CALL DSPLY(2) CALL GWflush(IR, -NOG) ENDIF * ENDDO CALL GWSETBK(IR, 1) CALL GWREFRESH(IR) CALL GWQUIT(IR) END C SUBROUTINE initial IMPLICIT REAL*8 (A-H,O-Z) REAL xmin, ymin, xmax, ymax, ch PARAMETER(KRED = 13, LSOLID = 1, NOTXORPEN = 11) COMMON A,B,T,ET,DT,XS,DX,ND,NPD,J COMMON /KDVPRM/AD12,BCD2,AD128,BCD COMMON /PRM/CR2I, CR2IM2 COMMON /PP/U(1000),WK1(1000),WK2(1000) COMMON /PLT/xmin, ymin, xmax, ymax, ch, NOG * CALL GWopen(IR,0) * CALL GWmode(IR, 2, 1) * CALL initwav * CR2I=1-DSQRT(.5D0) CR2IM2=2-CR2I AD12=DT*A*DX**(-1)/12 AD128=8*AD12 BCD2=DT*B*DX**(-3)/2 BCD=2*BCD2 * CALL wscale(XS,DX,U,ND) CALL GWanchor(IR, 1) CALL GWsetogn(NOG, 0, 0) NOG = NOG + 1 CALL GWsetogn(IR, 0, -NOG) CALL GWsetpen(IR, KRED, LSOLID, 0, NOTXORPEN) CALL gplot(XS,DX,U,ND) CALL GWflush(IR, -NOG) IF(MSGBOX('Ok ?') .NE. 1) THEN CALL GWquitx(IR,0) STOP ENDIF * END C SUBROUTINE initwav IMPLICIT REAL*8 (A-H,O-Z) DIMENSION NID(3), VD(2), WD(2), X0D(2) CHARACTER YN, str*80 LOGICAL DEFLT COMMON A,B,T,ET,DT,XS,DX,ND,NPD,J COMMON /II/TI,UI(1000),S1I,S2I,S3I COMMON /PP/U(1000),WK1(1000),WK2(1000) DATA V/-1D0/, NS/1/ * DATA XLD, NDD, XSD, DTD, ETD/40D0, 200, 0D0, 5D-3, 100D0/ DATA NID/2,2,0/ DATA VD/1D0,0.5D0/ DATA WD/2D0,2.828427125D0/ DATA X0D/10D0,25D0/ * CALL GWINPUT(IR,'use default parameters [Y/N] ? ', str) READ(str,'(A)') YN IF(YN .EQ. 'N' .OR. YN .EQ. 'n') THEN DEFLT = .FALSE. ELSE DEFLT = .TRUE. ENDIF * IF(DEFLT) THEN WRITE(*,'(/A)') 'K-dV Eq. : Ut+6*U*Ux+Uxxx=0' A=6 B=1 ELSE WRITE(*,'(/A)') 'K-dV Eq. : Ut+A*U*Ux+B*Uxxx=0' WRITE(*,'(/A,$)') + 'Enter parameters A[6.],B[1.] = ' READ(*,*) A,B ENDIF * IF(DEFLT) THEN XL=XLD ND=NDD XS=XSD WRITE(*,*) + 'width, # of div., origin of the frame =' + ,XL,ND,XS ELSE WRITE(*,*) + 'Enter width, # of div., org. of the frame =' READ (*,*) XL,ND,XS ENDIF * DX=XL/ND * IF(DEFLT) THEN DT=DTD T=0 ET=ETD WRITE(*,*) + 'time step, initial time, end time =' + ,DT,T,ET ELSE WRITE(*,'(/A,$)') + 'Enter DT(time slice),T(initial time),ET(end time) =' READ (*,*) DT,T,ET ENDIF * IF(DEFLT) THEN NPD=100 WRITE(*,*) + 'period for display =', NPD ELSE WRITE(*,'(/A,$)') + 'Enter period for display : NPD[100] =' READ (*,*) NPD IF(NPD.EQ.0) NPD=100 ENDIF * IF(DEFLT) THEN NI = NID(NS) ELSE WRITE(*,'(/A//A//A,$)') + ' 1. V*cos(PAI*x/W)', + ' 2. (V/2)*(sech(x/W))**2', + 'Select # of the initial wave form : ' READ(*,'(I1)') NI ENDIF DO WHILE(NI .EQ. 1 .OR. NI .EQ. 2) * IF(DEFLT) THEN V=VD(NS) W=WD(NS) X0=X0D(NS) WRITE(*,*) NS,'. parameters V,W,X0 = ' + , VD(NS),WD(NS),X0D(NS) ELSE IF(NI.EQ.1) THEN WRITE(*,'(/A)') 'INITIAL WAVE : V*cos(PAI*(x-X0)/W)' WRITE(*,'(/A,$)') 'Enter parameters V,W[1.],X0[0.] = ' READ(*,*) V,W,X0 IF(W.LE.0.) W=1 ELSE WRITE(*,'(/A)') 'INITIAL WAVE : (V/2)*(sech((x-X0)/W))**2' DO WHILE(V .LE. 0D0) WRITE(*,'(/A,$)') + 'Enter parameters V,W[2/|V|**(1/2)],X0 = ' READ(*,*) V,W,X0 ENDDO IF(W.LE.0.) W=2/DSQRT(DABS(V)) ENDIF ENDIF * WRITE(*,'(/A)') 'now making initial data' IF(NI.EQ.1) THEN DO I=1,ND U(I)=U(I)+V*DCOS(4*DATAN(1.D0)*(XS-X0+DX*(I-1))/W) ENDDO ELSE DO I=1,ND EXPX=DEXP(-2*DABS(XS-X0+DX*(I-1))/W) U(I)=U(I)+2*V*EXPX*(1+EXPX)**(-2) ENDDO ENDIF * NS = NS + 1 IF(DEFLT) THEN NI = NID(NS) ELSE WRITE(*,'(/A//A//A//A,$)') + ' 1. V*cos(PAI*x/W)', + ' 2. (V/2)*(sech(x/W))**2', + ' 3. End', + 'Select # of the initial wave form : ' READ(*,'(I1)') NI ENDIF ENDDO * XBS=XS TI=T J=0 DO I=1,ND UI(I)=U(I) WK1(I)=0 WK2(I)=U(I) ENDDO * CALL QCONST(U,S1I,S2I,S3I) S1=S1I S2=S2I S3=S3I WRITE(*,'(/3(A,F9.5))') + 'S1=',S1,' S2=',S2,' S3=',S3 * END C SUBROUTINE QCONST(U,S1,S2,S3) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION U(1000) COMMON A,B,T,ET,DT,XS,DX,ND,NPD,J S1=0 S2=0 S3=0 C=A/3/B D=1/DX/12 U2=U(ND-1) U3=U(ND) U4=U(1) U5=U(2) DO I=1,ND U1=U2 U2=U3 U3=U4 U4=U5 U5=U(MOD(I+1,ND)+1) CS=(1+MOD(I,2))*DX/3 S1=S1+CS*U3 S2=S2+CS*U3*U3 S3=S3+CS*(C*U3**3-(D*(U1-U5+8*(U4-U2)))**2) ENDDO END C SUBROUTINE DSPLY(NW) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER buff*60 REAL xmin, ymin, xmax, ymax, ch PARAMETER(KBLUE = 16, LSOLID = 1, NOTXORPEN = 11) COMMON A,B,T,ET,DT,XS,DX,ND,NPD,J COMMON /II/TI,UI(1000),S1I,S2I,S3I COMMON /PP/U(1000),WK1(1000),WK2(1000) COMMON /PLT/xmin, ymin, xmax, ymax, ch, NOG * WRITE(buff,'(A,F8.4,A,I6,A)') 't =',T,' ( step =',J,' )' CALL GWputtxt(IR, xmin, ymax*1.01, buff) * CALL QCONST(U,S1,S2,S3) WRITE(buff,'(1P,3(A,E10.2))') + 'err1 =',S1-S1I,', err2 =',S2-S2I,', err3 =',S3-S3I CALL GWputtxt(IR, xmin, ymin-ch*1.1, buff) CALL GWsetpen(IR, KBLUE, LSOLID, 0, NOTXORPEN) IF(NW .EQ. 2) THEN CALL gplot(XS,DX,U,ND) TI = T * S1I = S1 * S2I = S2 * S3I = S3 ELSE CALL gplot(XS,DX,U,ND) ENDIF END C SUBROUTINE gplot( XS,DX,U,ND) REAL*8 XS,DX,U(ND) REAL POINTS(2,0:1000) REAL xmin, ymin, xmax, ymax, ch COMMON /PLT/xmin, ymin, xmax, ymax, ch, NOG * POINTS(1,0)=XS POINTS(2,0)=U(ND) DO id = 1,ND POINTS(1,id)=XS+DX*id POINTS(2,id)=U(id) ENDDO CALL GWpolylin(IR, POINTS, ND+1) * END C SUBROUTINE wscale(XS, DX, U, ND) REAL*8 XS, DX, U(ND) REAL xmin, ymin, xmax, ymax, ch PARAMETER(KBLACK = 0, LSOLID = 1, LDASH = 2, MCOPYPEN = 4) COMMON /PLT/xmin, ymin, xmax, ymax, ch, NOG * width = DX*ND xmin = XS xmax = XS+width ymin = U(1) ymax = U(1) do id = 1,ND if( U(id) .lt. ymin ) ymin = U(id) if( U(id) .gt. ymax ) ymax = U(id) ENDDO * height = ymax - ymin ymax = ymax + height/4 ymin = ymin - height/4 ch = height/20 * CALL GWclear(IR,-1) CALL GWindow(IR, xmin-width*0.1,ymin-height*0.1, + xmax+width*0.1,ymax+height*0.1 ) CALL GWsettxt(IR, ch, 0.0, -1, -1, -1, ' ') * CALL GWsetpen(IR,KBLACK, LSOLID, 0, MCOPYPEN) * CALL GWline(IR, xmin,ymin,xmax,ymin) CALL GWline(IR, xmin,ymax,xmax,ymax) CALL GWline(IR, xmin,ymin,xmin,ymax) CALL GWline(IR, xmax,ymin,xmax,ymax) * CALL GWsetpen(IR, KBLACK, LDASH, 0, MCOPYPEN) * pdy = (ymax-ymin)/10. DO n = 1,9 CALL GWline(IR, xmin,ymin+n*pdy,xmax,ymin+n*pdy) ENDDO * pdx = (xmax-xmin)/10. DO n = 1,9 CALL GWline(IR, xmin+n*pdx,ymin,xmin+n*pdx,ymax) ENDDO * END C SUBROUTINE UPDATE IMPLICIT REAL*8 (A-H,O-Z) DIMENSION WK3(1000) COMMON A,B,T,ET,DT,XS,DX,ND,NPD,J COMMON /KDVPRM/AD12,BCD2,AD128,BCD COMMON /PP/U(1000),WK1(1000),WK2(1000) COMMON /PRM/CR2I, CR2IM2 U2=WK2(ND-1) U3=WK2(ND) U4=WK2(1) U5=WK2(2) DO K=1,ND U1=U2 U2=U3 U3=U4 U4=U5 U5=WK2(MOD(K+1,ND)+1) T2= (BCD2 - AD12*U3)*(U1-U5) + (BCD - AD128*U3)*(U4-U2) T3=T2/2-WK1(K) WK3(K)=WK2(K)+T3 T3=WK3(K)-WK2(K) WK1(K)=WK1(K)+3*T3-T2/2 ENDDO * U2=WK3(ND-1) U3=WK3(ND) U4=WK3(1) U5=WK3(2) DO K=1,ND U1=U2 U2=U3 U3=U4 U4=U5 U5=WK3(MOD(K+1,ND)+1) T2= (BCD2 - AD12*U3)*(U1-U5) + (BCD - AD128*U3)*(U4-U2) T3=CR2I*(T2-WK1(K)) WK2(K)=WK3(K)+T3 T3=WK2(K)-WK3(K) WK1(K)=WK1(K)+3*T3-CR2I*T2 ENDDO * U2=WK2(ND-1) U3=WK2(ND) U4=WK2(1) U5=WK2(2) DO K=1,ND U1=U2 U2=U3 U3=U4 U4=U5 U5=WK2(MOD(K+1,ND)+1) T2= (BCD2 - AD12*U3)*(U1-U5) + (BCD - AD128*U3)*(U4-U2) T3=CR2IM2*(T2-WK1(K)) WK3(K)=WK2(K)+T3 T3=WK3(K)-WK2(K) WK1(K)=WK1(K)+3*T3-CR2IM2*T2 ENDDO * U2=WK3(ND-1) U3=WK3(ND) U4=WK3(1) U5=WK3(2) DO K=1,ND U1=U2 U2=U3 U3=U4 U4=U5 U5=WK3(MOD(K+1,ND)+1) T2= (BCD2 - AD12*U3)*(U1-U5) + (BCD - AD128*U3)*(U4-U2) T3=(T2-2*WK1(K))/6 WK2(K)=WK3(K)+T3 U(K)=WK2(K) T3=WK2(K)-WK3(K) WK1(K)=WK1(K)+3*T3-T2/2 ENDDO * J=J+1 T=T+DT END