c c Sample_Code_11 c Solves a linear Parabolic PDE using Chebyshev collocation Method c c dy/dt = d2y/dx2 c t = 0 y = 1 c x = 0 dy/dx = 0 c x = 1 y = 1 c c IMETHD = 1 Exp. Euler c IMETHD = 2 Imp. Euler c IMPLICIT DOUBLE PRECISION (a-h,o-z) parameter (N=16,N1=N+1) DOUBLE PRECISION YOLD(N1) DOUBLE PRECISION A(N1,N1),B(N1) integer ipvt(N1) COMMON /CHEB/ x1,x2,PI,Alpha,Beta COMMON /TRANS/ C(N1),CBAR(N1),CT(N1,N1),CTINV(N1,N1),X(N1) COMMON /COLLEC/ D(N1,N1),DD(N1,N1) COMMON /PARAM/ Phi,Bi COMMON /PARAM2/ TMAX,H,IMETHD,ISTEP OPEN(11,FILE='sample_code_11_output1.txt') OPEN(12,FILE='sample_code_11_output2.txt') c IMETHD = 1 PI = 4.d0*DATAN(1.d0) x1 = 0.d0 x2 = 1.d0 Phi = 1.d0 Bi = 1.d0 h = 1.d-4 TMAX = 10.d0 IPRINT = (TMAX/h)/100 c call Chebyshev c c Provide Initial Solution c do i = 1, N1 YOLD(i) = 1.d0 enddo c c Print the Initial solution write(11,*) write(11,'(3x,A,f10.6)') 'Time =',Time do i = 1, N1 write(11,'(10x,2f16.6)') X(i),YOLD(i) enddo c c Evaluate Matrix A once (linear equations) call MATRIX(A) c c LU-decompose matrix A c Note that matrix A will contains on return the L and U matrixes (original is destroyed) call ludcmp(A,N1,N1,ipvt,dummy) c c Start Integration Loop count = 0.d0 istep = 0 Time = 0.d0 10 Time = Time + h istep = istep + 1 c c Set Termination Criteria for the Integration Loop if((Time-TMAX) .gt. 1.d-14) stop c c Evaluate RHS of linear equations call FUNC(YOLD,B) c c Solve the system of linear equations call lubksb(A,N1,N1,ipvt,B) c c Print the solution for the current time if (istep .eq. iprint) then write(11,'(3x,A,f10.6)') 'Time =',Time do i = 1, N1 write(11,'(10x,2f16.6)') X(i),B(i) enddo write(12,'(2f16.6)') Time,B(N1-1) write( *,'(2f16.6)') Time,B(N1-1) istep = 0 endif c do i = 1, N1 YOLD(i) = B(i) enddo c go to 10 c 1000 stop end c c c c c c c c c c c c SUBROUTINE MATRIX(A) IMPLICIT DOUBLE PRECISION (a-h,o-z) parameter (N=16,N1=N+1) DOUBLE PRECISION A(N1,N1) COMMON /CHEB/ x1,x2,PI,Alpha,Beta COMMON /TRANS/ C(N1),CBAR(N1),CT(N1,N1),CTINV(N1,N1),X(N1) COMMON /COLLEC/ D(N1,N1),DD(N1,N1) COMMON /PARAM2/ TMAX,H,IMETHD,ISTEP c do i = 1, N1 do j = 1, N1 A(i,j) = 0.d0 enddo enddo c c Construct Matrix using Chebyshev Collocation c If (IMETHD .eq. 1) then c c Explicit Euler Integration do i = 2, N1-1 A(i,i) = 1.d0 enddo go to 100 endif If (IMETHD .eq. 2) then c c Implicit Euler Integration do i = 2, N1-1 A(i,i) = 1.d0 do j = 1, N1 A(i,j) = A(i,j)-h*DD(i,j) enddo enddo go to 100 endif c c Implement boundary conditions c 100 do j = 1, N1 A(1,j) = D(1,j) enddo A(N1,N1) = 1.d0 c return end c c c c c c c c c c c c c Subroutine Func to Input RHS functions of PDE SUBROUTINE FUNC(YOLD,B) IMPLICIT DOUBLE PRECISION (a-h,o-z) parameter (N=16,N1=N+1) DOUBLE PRECISION YOLD(N1),B(N1) COMMON /CHEB/ x1,x2,PI,Alpha,Beta COMMON /TRANS/ C(N1),CBAR(N1),CT(N1,N1),CTINV(N1,N1),X(N1) COMMON /COLLEC/ D(N1,N1),DD(N1,N1) COMMON /PARAM/ Phi,Bi COMMON /PARAM2/ TMAX,H,IMETHD,ISTEP c do i = 1, N1 B(i) = 0.d0 enddo c c Construct Residuals for the RHS functions of PDE using Chebyshev Collocation c If (IMETHD .eq. 1) then c Explicit Euler Integration do i = 2, N1-1 B(i) = YOLD(i) do j = 1, N1 B(i) = B(i) + h*DD(i,j)*YOLD(j) enddo enddo go to 100 endif c If (IMETHD .eq. 2) then c Implicit Euler Integration do i = 2, N1-1 B(i) = YOLD(i) enddo go to 100 endif c c Boundary Conditions c In Collecation, the first equation for the first BC and the last equation for the second BC 100 B( 1) = 0.d0 B(N1) = 0.d0 c return end c c c SUBROUTINE Chebyshev IMPLICIT DOUBLE PRECISION (a-h,o-z) parameter (N=16,N1=N+1) DOUBLE PRECISION XCHEB(N1) COMMON /CHEB/ x1,x2,PI,Alpha,Beta COMMON /TRANS/ C(N1),CBAR(N1),CT(N1,N1),CTINV(N1,N1),X(N1) COMMON /COLLEC/ D(N1,N1),DD(N1,N1) c c Vector XCHEB (1,-1) contains the Gauss-Lobatto grid. do i = 1, N1 XCHEB(i) = dcos((i-1)*pi/N) CBAR(i) = 1.d0 C(i) = 1.d0 enddo c C(1) = 2.d0 CBAR(1) = 2.d0 CBAR(N1) = 2.d0 c c Vector X maps XCHEB to the interval (x1,x2) ==> X(i) = Alpha + Beta*XCHEB(i) Alpha = (x1+x2)/2.d0 Beta = (x1-x2)/2.d0 c do i = 1, N1 X(i) = Alpha + Beta*XCHEB(i) enddo c do i = 1, N1 do j = 1, N1 CT(i,j) = 2.d0/N/CBAR(i)/CBAR(j)*dcos(pi*(i-1)*(j-1)/N) CTINV(i,j) = dcos(pi*(i-1)*(j-1)/N) enddo enddo c do i = 1, N1 ii = i-1 do j = 1, N1 jj = j-1 if (ii .ne. jj) then D(i,j) = CBAR(i)/CBAR(j)*(-1)**(ii+jj)/(XCHEB(i)-XCHEB(j)) else if (ii .ne. 0 .and. ii .ne. N) then D(i,j) = -XCHEB(j)/2.d0/(1.d0-XCHEB(j)**2) endif endif enddo enddo D( 1, 1) = (2.d0*dble(N)**2+1.d0)/6.d0 D(N1,N1) =-(2.d0*dble(N)**2+1.d0)/6.d0 c do i=1,N1 do j=1,N1 DD(i,j) = 0.d0 do k=1,N1 DD(i,j) = DD(i,j) + D(i,k)*D(k,j) enddo enddo enddo c do i = 1, N1 do j = 1, N1 D(i,j) = D(i,j)/BETA DD(i,j) = DD(i,j)/BETA/BETA enddo enddo c return end c c c c c c c c c c c c c c SUBROUTINE ludcmp(a,n,np,indx,d) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER n,np,indx(n),NMAX DOUBLE PRECISION d,a(np,np),TINY PARAMETER (NMAX=500,TINY=1.0e-20) INTEGER i,imax,j,k DOUBLE PRECISION aamax,dum,sum,vv(NMAX) d=1. do 12 i=1,n aamax=0.d0 do 11 j=1,n if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j)) 11 continue if (aamax.eq.0.d0) pause 'singular matrix in ludcmp' vv(i)=1./aamax 12 continue do 19 j=1,n do 14 i=1,j-1 sum=a(i,j) do 13 k=1,i-1 sum=sum-a(i,k)*a(k,j) 13 continue a(i,j)=sum 14 continue aamax=0.d0 do 16 i=j,n sum=a(i,j) do 15 k=1,j-1 sum=sum-a(i,k)*a(k,j) 15 continue a(i,j)=sum dum=vv(i)*abs(sum) if (dum.ge.aamax) then imax=i aamax=dum endif 16 continue if (j.ne.imax)then do 17 k=1,n dum=a(imax,k) a(imax,k)=a(j,k) a(j,k)=dum 17 continue d=-d vv(imax)=vv(j) endif indx(j)=imax if(a(j,j).eq.0.d0)a(j,j)=TINY if(j.ne.n)then dum=1./a(j,j) do 18 i=j+1,n a(i,j)=a(i,j)*dum 18 continue endif 19 continue return END c SUBROUTINE lubksb(a,n,np,indx,b) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER n,np,indx(n) DOUBLE PRECISION a(np,np),b(n) INTEGER i,ii,j,ll DOUBLE PRECISION sum ii=0 do 12 i=1,n ll=indx(i) sum=b(ll) b(ll)=b(i) if (ii.ne.0)then do 11 j=ii,i-1 sum=sum-a(i,j)*b(j) 11 continue else if (sum.ne.0.d0) then ii=i endif b(i)=sum 12 continue do 14 i=n,1,-1 sum=b(i) do 13 j=i+1,n sum=sum-a(i,j)*b(j) 13 continue b(i)=sum/a(i,i) 14 continue return END