c c Sample_Code_6 c Solves a BVP using Chebyshev Galerkin Method c IMPLICIT DOUBLE PRECISION (a-h,o-z) parameter (N=10,N1=N+1) DOUBLE PRECISION A(N1,N1),B(N1),Y(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) OPEN(11,FILE='sample_code_6_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,d) c do i = 1, N1-2 B(i) = 0.d0 enddo B(N1-1) = 0.d0 B(N1) = 1.d0 c call lubksb(A,N1,N1,ipvt,B) call SPRE(B,Y) 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),Y(i),ye write(* ,'(3f16.6)') Z(i),Y(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) c do i = 1, N1 do j = 1, N1 A(i,j) = 0.d0 enddo enddo c c Construct the residals in the spectral domain c do i = 1, N1-2 ri = dble(i) - 1.d0 c A(i,i) = A(i,i) + 4.d0 c 1st-oder derivatives do ip = i + 1, N1, 2 rp = dble(ip) - 1.d0 A(i,ip) = A(i,ip) + 5.d0*(2.d0/C(i)*rp)/Beta enddo c 2nd-oder derivatives do ip = i + 2, N1, 2 rp = dble(ip) - 1.d0 A(i,ip) = A(i,ip) + 1.d0/C(i)*rp*(rp*rp-ri*ri)/Beta/Beta enddo enddo c c Implement boundary conditions c do j = 1, N1 A(N1-1,j) = 1.d0 A(N1 ,j) = (-1.d0)**(j-1) enddo 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) 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 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