c c Sample_Code_7 c Solves a nonlinear BVP using Chebyshev Galerkin Method c IMPLICIT DOUBLE PRECISION (a-h,o-z) parameter (N=16,N1=N+1) DOUBLE PRECISION Y(N1),YHAT(N1),YHAT_G(N1) COMMON /CHEB/ z1,z2,PI,Alpha,Beta COMMON /TRANS/ C(N1),CBAR(N1),CT(N1,N1),CTINV(N1,N1),Z(N1) COMMON /PARAM/ Phi,Bi OPEN(11,FILE='sample_code_7_output.txt') c PI = 4.d0*DATAN(1.d0) z1 = 0.d0 z2 = 1.d0 Phi = 1.d0 Bi = 1.d0 c call Chebyshev c do i = 1, N1 Y(i) = 1.d0/(dexp(-1.d0)-dexp(-4.d0))* ! (dexp(-Z(i))-dexp(-4.d0*Z(i)))*0.5 enddo c call resp(Y,YHAT_G) c call Newton(YHAT_G,YHAT) c call spre(YHAT,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,'(3f20.14)') Z(i),Y(i),ye write(* ,'(3f20.14)') Z(i),Y(i),ye enddo c 1000 stop end c c c c c c c c Subroutine Newton to to solve a system of algebraic equations iteratively subroutine Newton(UNKNON_G,UNKNON) IMPLICIT DOUBLE PRECISION (A-H,O-Z) parameter (N=16,N1=N+1) DOUBLE PRECISION A(N1,N1),UNKNON_G(N1),UNKNON(N1),B(N1) integer ipvt(N1) c itmax = 10 tol = 1.d-12 c c initialize matrix A and vector B do i = 1, N1 UNKNON(i) = UNKNON_G(i) B(i) = 0.d0 do j = 1, N1 A(i,j)=0.d0 enddo enddo c c Start Iterative Solution iter = 0 10 iter = iter + 1 if (iter .gt. itmax) go to 12 c c Input system of nonlinear equations c on return B contains the residuals call residuals(UNKNON,B) c c Check the norm of the residuals rnorm = 0.d0 do i = 1, N1 rnorm = rnorm + B(i)**2 enddo rnorm = dsqrt(rnorm) write(*,'(3x,a,i3,5x,a,e17.10)') ! 'iteration #',iter, 'Tolerance =', rnorm if (rnorm .lt. tol) return c c Evaluate the Jacobian Matrix c IFLAG = 0 ==> analatical jacobian c IFLAG = 1 ==> numerical jacobian IFLAG = 1 call jacobian(IFLAG,UNKNON,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 Solve the equations: A .(Xiter+1-Xiter) =-B, by backsubstitution c Note that vector B will contains the unknowns X (original is destroyed) c for Newton Raphson method B = Xiter+1 - Xiter do i = 1, N1 B(i) = -B(i) enddo c call lubksb(A,N1,N1,ipvt,B) c do i= 1, N1 UNKNON(i) = UNKNON(i)+B(i) enddo go to 10 12 write(11,*) 'No Convergence !' write(*,*) 'No Convergence !' stop return end c c c c c c c c c c c c c Subroutine Residuals to Input System of Non-Linear Equations SUBROUTINE RESIDUALS(YHAT,RHAT) IMPLICIT DOUBLE PRECISION (a-h,o-z) parameter (N=16,N1=N+1) DOUBLE PRECISION Y(N1),YHAT(N1),R(N1),RHAT(N1) DOUBLE PRECISION D1YHAT(N1),D2YHAT(N1),D1Y(N1),D2Y(N1) COMMON /CHEB/ z1,z2,PI,Alpha,Beta COMMON /TRANS/ C(N1),CBAR(N1),CT(N1,N1),CTINV(N1,N1),Z(N1) COMMON /PARAM/ Phi,Bi c do i = 1, N1 R(i) = 0.d0 RHAT(i) = 0.d0 enddo c c Evaluate the derivatines in the spectral domain c call Derivative(YHAT,D1YHAT,D2YHAT) c call spre(YHAT,Y) call spre(D1YHAT,D1Y) call spre(D2YHAT,D2Y) c c Construct Residuals in the Physical Domain then Transform to Spectral Domain c do i = 1, N1 R(i) = (D2Y(i)+5.d0*D1Y(i)+4.d0*Y(i))*Y(i) enddo c c Boundary Conditions c In physical domain, the first equation for the first BC and the last equation for the second BC c R(1 ) = Y(1 )-0.d0 R(N1) = Y(N1)-1.d0 c c Transform the residuals back to the spectral domain c call resp(R,RHAT) c return end c c c c c c c c c Subroutine jacobian to Evaluate the Jacobian Matrix subroutine jacobian(IFLAG,UNKNON,A) IMPLICIT DOUBLE PRECISION (A-H,O-Z) parameter (N=16,N1=N+1) double precision A(N1,N1),UNKNON(N1),UNKNONJ(N1),FUN(N1),FUNJ(N1) parameter (epsl=1.d-5) c do i = 1, N1 do j = 1, N1 A(i,j)=0.d0 enddo enddo c if (IFLAG .gt. 0) go to 10 c c Evaluate the jacobian matrix analatically return c c Evaluate the jacobian matrix numerically 10 do i = 1, N1 UNKNONJ(i)=UNKNON(i) enddo call residuals(UNKNON,FUN) do i = 1, N1 diff = dmax1(epsl,dabs(epsl*UNKNON(i))) UNKNONJ(i) = UNKNON(i)+diff call residuals(UNKNONJ,FUNJ) do j = 1, N1 A(j,i) = (FUNJ(j)-FUN(j))/diff enddo UNKNONJ(i)=UNKNON(i) enddo 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/ 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 c c c SUBROUTINE Derivative(YHAT,D1YHAT,D2YHAT) IMPLICIT DOUBLE PRECISION (a-h,o-z) parameter (N=16,N1=N+1) DOUBLE PRECISION YHAT(N1),D1YHAT(N1),D2YHAT(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 D1YHAT(i) = 0.d0 D2YHAT(i) = 0.d0 enddo c do i = 1, N1 ri = dble(i) - 1.d0 c 1st-oder derivatives do ip = i + 1, N1, 2 rp = dble(ip) - 1.d0 D1YHAT(i) = D1YHAT(i) + (2.d0/C(i)*rp)/Beta*YHAT(ip) enddo c 2nd-oder derivatives do ip = i + 2, N1, 2 rp = dble(ip) - 1.d0 D2YHAT(i) = D2YHAT(i) + ! 1.d0/C(i)*rp*(rp*rp-ri*ri)/Beta/Beta*YHAT(ip) enddo enddo c return end cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx subroutine RESP(AR,AS) cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx implicit double precision (a-h,o-z) parameter (N=16,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=16,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