c c Sample_Code_8 c Solves a linear BVP using Chebyshev Collocation Method c IMPLICIT DOUBLE PRECISION (a-h,o-z) parameter (N=10,N1=N+1) DOUBLE PRECISION A(N1,N1),B(N1) INTEGER IPVT(N1) COMMON /CHEB/ z1,z2,PI,Alpha,Beta COMMON /TRANS/ C(N1),CBAR(N1),CT(N1,N1),CTINV(N1,N1),Z(N1) COMMON /COLLEC/ D(N1,N1),DD(N1,N1) OPEN(11,FILE='sample_code_8_output.txt') c PI = 4.d0*ATAN(1.d0) z1 = 0.d0 z2 = 1.d0 c call Chebyshev c call MATRIX(A) c call ludcmp(A,N1,N1,ipvt,dummy) c do i = 1, N1 B(i) = 0.d0 enddo B(N1) = 1.d0 c call lubksb(A,N1,N1,ipvt,B) c do i = 1, N1 ye = 1.d0/(dexp(-1.d0)-dexp(-4.d0))* ! (dexp(-Z(i))-dexp(-4.d0*Z(i))) write(11,'(3f16.6)') Z(i),B(i),ye write(* ,'(3f16.6)') Z(i),B(i),ye enddo c 1000 stop end c c c SUBROUTINE MATRIX(A) IMPLICIT DOUBLE PRECISION (a-h,o-z) parameter (N=10,N1=N+1) DOUBLE PRECISION A(N1,N1) COMMON /CHEB/ z1,z2,PI,Alpha,Beta COMMON /TRANS/ C(N1),CBAR(N1),CT(N1,N1),CTINV(N1,N1),Z(N1) COMMON /COLLEC/ D(N1,N1),DD(N1,N1) c do i = 1, N1 do j = 1, N1 A(i,j) = 0.d0 enddo enddo c c Construct Residuals using Chebyshev Collocation c do i = 2, N1-1 A(i,i) = 4.d0 do j = 1, N1 A(i,j) = A(i,j)+DD(i,j)+5.d0*D(i,j) enddo enddo c c Implement boundary conditions c A( 1, 1) = 1.d0 A(N1,N1) = 1.d0 c return end c c c SUBROUTINE Chebyshev IMPLICIT DOUBLE PRECISION (a-h,o-z) parameter (N=10,N1=N+1) DOUBLE PRECISION XCHEB(N1) COMMON /CHEB/ z1,z2,PI,Alpha,Beta COMMON /TRANS/ C(N1),CBAR(N1),CT(N1,N1),CTINV(N1,N1),Z(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 (z1,z2) ==> X(i) = Alpha + Beta*XCHEB(i) Alpha = (z1+z2)/2.d0 Beta = (z1-z2)/2.d0 c do i = 1, N1 Z(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 cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx subroutine RESP(AR,AS) cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx implicit double precision (a-h,o-z) parameter (N=10,N1=N+1) DOUBLE PRECISION AR(N1),AS(N1) COMMON /TRANS/ C(N1),CBAR(N1),CT(N1,N1),CTINV(N1,N1),Z(N1) c c invert chebychev c do i = 1, N1 AS(i) = 0.d0 do j = 1, N1 AS(i) = AS(i) + CT(i,j) * AR(j) enddo enddo c return end cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx subroutine SPRE(AS,AR) cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx implicit double precision (a-h,o-z) parameter (N=10,N1=N+1) DOUBLE PRECISION AR(N1),AS(N1) COMMON /TRANS/ C(N1),CBAR(N1),CT(N1,N1),CTINV(N1,N1),Z(N1) c c invert chebychev c do i = 1, N1 AR(i) = 0.d0 do j = 1, N1 AR(i) = AR(i) + CTINV(i,j) * AS(j) 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