parameter (N=64, NN=2*N+1) implicit real*8(a-h,o-z) real x(0:NN+1),ubar(0:NN+1,0:3) real uapp1(0:NN+1), xL(0:NN+1,0:3) complex abar(0:NN+2), al(0:NN+2), si(0:NN+1) complex w, A, y, ss, ss1, ss2, sum, sum1 pi=4.0*atan(1.0) w=cmplx(0, 1) alpha=15.0*log(10.0) dx=2.0*pi/NN dx1=0.5*dx dt=0.001 ttime=0.0 Tmax=0.2 par=1.9 c compute the initial fourier coefficient do j=0, NN x(j)=j*dx if ( (x(j)-dx1) .lt. par .and. (x(j)+dx1) .gt. par) jp=j enddo do j=0, NN x(j)=j*dx c ubar(j, 0)=0.3+0.7*(-cos(x(j)+0.5*dx)+cos(x(j)-0.5*dx))/dx s1=cos(0.5*(x(j)+dx1)) s2=cos(0.5*(x(j)-dx1)) if (j .lt. jp) then ubar(j, 0)=-2.0*(s1-s2)/dx elseif (j .gt. jp) then ubar(j, 0)=2.0*(s1-s2)/dx else p1=-2.0*(cos(0.5*par)-s2)/dx p2=2.0*(s1-cos(0.5*par))/dx ubar(j, 0)=p1+p2 endif enddo 99 do NR=0, 2 do k=0, NN-1 k1=k-0.5*(NN-1) ss=0.0 do j=0, NN-1 s1=cos(k1*x(j)) s2=-sin(k1*x(j)) ss=ss+ubar(j,NR)*cmplx(s1,s2) enddo abar(k)=ss/NN enddo do k=0, NN-1 k1=k-0.5*(NN-1) if (k1 .eq. 0) then si(k)=1.0 else si(k)=sin(0.5*k1*dx)/(0.5*k1*dx) endif al(k)=abar(k)/si(k) enddo do j=0, NN-1 xj=(j+0.5)*dx sum=0.0 do k=1, N tau=1.0 ps=k*xj s1=cos(ps) s2=sin(ps) ss2=cmplx(s1, s2) ss=al(k+N)*ss2 sum=sum+tau*ss k1=-k ps=k1*xj s1=cos(ps) s2=sin(ps) ss2=cmplx(s1, s2) ss=al(k1+N)*ss2 sum=sum+tau*ss enddo uapp1(j)=real(al(N)+sum) enddo do j=0, NN-1 ss1=uapp1(j) if (j .eq. 0) then ss2=uapp1(NN-1) else ss2=uapp1(j-1) endif xL(j, NR)=-0.5*(ss1*ss1-ss2*ss2)/dx if (NR .eq. 0) then ubar(j,1)=ubar(j,0)+dt*xL(j,0) endif if (NR .eq. 1) then sa=0.75*ubar(j,0)+0.25*ubar(j,1)+0.25*dt*xL(j,1) ubar(j,2)=sa endif if (NR .eq. 2) then sa=ubar(j,0)/3.0+ubar(j,2)*2/3.0+dt*xL(j,2)*2/3.0 ubar(j,3)=sa endif enddo enddo ttime=ttime+dt print *, ttime, ' time' if (ttime .lt. Tmax) then do j=0, NN ubar(j,0)=ubar(j, 3) enddo goto 99 endif do k=0, NN-1 k1=k-0.5*(NN-1) ss=0.0 do j=0, NN-1 p1=k1*x(j) s1=cos(p1) s2=-sin(p1) ss=ss+ubar(j,NR)*cmplx(s1,s2) enddo abar(k)=ss/NN al(k)=abar(k)/si(k) enddo do j=0, NN-1 xj=x(j) sum=0.0 do k=1, N tau=1.0 ps=k*xj s1=cos(ps) s2=sin(ps) ss2=cmplx(s1, s2) ss=al(k+N)*ss2 sum=sum+tau*ss k1=-k ps=k1*xj s1=cos(ps) s2=sin(ps) ss2=cmplx(s1, s2) ss=al(k1+N)*ss2 sum=sum+tau*ss enddo uapp1(j)=real(al(N)+sum) enddo write (17, *) 'x=[' do j=0, NN-1 write (17, *) x(j) enddo write (17, *) '];' write (18, *) 'y=[' do j=0, NN-1 write (18, *) uapp1(j) enddo write (18, *) '];' print *, A, y, ' A, y' stop end