c c Sample_Code_4 c Solves a single BVP using the shooting method c IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (N = 2) DOUBLE PRECISION Y(N),YOLD(N) COMMON x1,x2,h,phi1,phi2 OPEN(11,FILE='sample_code_4_output.txt') c tol = 1.d-10 ! Maximum Tolerance itmax = 20 ! Maximum Number of Iterations IFLAG = 1 ! Flag for the Jacobian Matrix (0 -> analytical Jacobian, 1 -> Numerical Jacobian) x1 = 0.d0 ! Starting Time x2 = 1.d0 ! Ending Time h = 5.d-2 ! Time Step phi1 = 1.d0 ! Parameter phi2 = 5.d-1 ! Parameter IFLAG2 = 4 ! 1-> Adams Bashforth, 2-> Predictor Corrector, 3-> Adams Moulten, 4-> RKG4 epsl = 1.d-5 ! Epsilon (a small number) c c provide initial condition for YOLD(N) Beta = 2.d0 znm1 = 0.1d0 zn = znm1+epsl c c Integrate IVP's for the first guess of IC YOLD(1) = 1.d0 YOLD(2) = znm1 c IF (IFLAG2 .eq. 1) GO TO 10 IF (IFLAG2 .eq. 2) GO TO 20 IF (IFLAG2 .eq. 3) GO TO 30 IF (IFLAG2 .eq. 4) GO TO 40 c c Integrate IVP's using 2nd-order Adams-Bashforth Method 10 call Adams_Bashforth(N,YOLD,Y) go to 50 c c Integrate IVP's using 2nd-order Predictor Corrector Method 20 call Predictor_Corrector(N,YOLD,Y) go to 50 c c Integrate IVP's using 2nd-order Adams-Moulten Method 30 call Adams_Moulten(N,ITMAX,TOL,IFLAG,YOLD,Y) go to 50 c c Integrate IVP's using 4th order Runge-Kutta-Gill Method 40 call RKG4(N,YOLD,Y) go to 50 c 50 phinm1 = Y(1) c c Integrate IVP's for the second guess of IC YOLD(1) = 1.d0 YOLD(2) = zn c IF (IFLAG2 .eq. 1) GO TO 11 IF (IFLAG2 .eq. 2) GO TO 21 IF (IFLAG2 .eq. 3) GO TO 31 IF (IFLAG2 .eq. 4) GO TO 41 c c Integrate IVP's using 2nd-order Adams-Bashforth Method 11 call Adams_Bashforth(N,YOLD,Y) go to 51 c c Integrate IVP's using 2nd-order Predictor Corrector Method 21 call Predictor_Corrector(N,YOLD,Y) go to 51 c c Integrate IVP's using 2nd-order Adams-Moulten Method 31 call Adams_Moulten(N,ITMAX,TOL,IFLAG,YOLD,Y) go to 51 c c Integrate IVP's using 4th order Runge-Kutta-Gill Method 41 call RKG4(N,YOLD,Y) go to 51 c 51 phin = Y(1) c c c c c c c Shooting Method Loop c c c c c c ishoot = 0 60 ishoot = ishoot + 1 if (ishoot .gt. itmax) stop write(*,'(5x,a,i3,5x,a,f20.12)') ! 'iteration =',ishoot,'error =',phin-beta c c Set Termination Criteria for the Shooting Method Loop if (dabs(phin-beta) .lt. tol) stop znp1 = zn + (beta-phin)*(zn-znm1)/(phin-phinm1) c c Integrate IVP's for the current guess of IC YOLD(1) = 1.d0 YOLD(2) = znp1 c IF (IFLAG2 .eq. 1) GO TO 12 IF (IFLAG2 .eq. 2) GO TO 22 IF (IFLAG2 .eq. 3) GO TO 32 IF (IFLAG2 .eq. 4) GO TO 42 c c Integrate IVP's using 2nd-order Adams-Bashforth Method 12 call Adams_Bashforth(N,YOLD,Y) go to 52 c c Integrate IVP's using 2nd-order Predictor Corrector Method 22 call Predictor_Corrector(N,YOLD,Y) go to 52 c c Integrate IVP's using 2nd-order Adams-Moulten Method 32 call Adams_Moulten(N,ITMAX,TOL,IFLAG,YOLD,Y) go to 52 c c Integrate IVP's using 4th order Runge-Kutta-Gill Method 42 call RKG4(N,YOLD,Y) go to 52 52 znm1 = zn zn = znp1 phinm1 = phin phin = Y(1) c go to 60 c end c c c c Subroutine Func to Input RHS functions of IVP c c c c subroutine Func(N,Y,FUN) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION Y(N),FUN(N) COMMON x1,x2,h,phi1,phi2 c do i = 1, N FUN(i) = 0.d0 enddo c FUN(1) = Y(2) FUN(2) = Y(1) c return end c c c c Subroutine Adams_Bashforth to integrate IVP's c c c subroutine Adams_Bashforth(N,YOLD,Y) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION Y(N),YOLD(N),FUNOLD1(N),FUNOLD2(N) COMMON x1,x2,h,phi1,phi2 c istep = 0 Time = x1 c c Print the Initial solution to output file write(11,*) write(11,'(12x,a,13x,a,13x,a)') 'X','Y(1)','Y(2)' write(11,'(3f16.6)') Time,YOLD(1),YOLD(2) c c Start Integration Loop 10 Time = Time + h istep = istep + 1 c c Set Termination Criteria for the Integration Loop if((Time-x2) .gt. 1.d-14) return c Evaluate the RHS Function for previous time call Func(N,YOLD,FUNOLD1) c c For the first time step apply 1st-order explicit Euler if (istep .eq. 1) then do i = 1, N Y(i) = YOLD(i)+h*FUNOLD1(i) enddo else c c Apply 2nd-order AB Formula for Remaining Time Steps do i = 1, N Y(i) = YOLD(i)+h/2.d0*(3.d0*FUNOLD1(i)-FUNOLD2(i)) enddo endif c c Print the solution for the current time step to output file write(11,'(3f16.6)') Time,Y(1),Y(2) c c Store Previous Time Values do i = 1, N YOLD(i) = Y(i) FUNOLD2(i) = FUNOLD1(i) enddo c go to 10 c end c c c c Subroutine Predictor-Corrector to integrate IVP's c c c subroutine Predictor_Corrector(N,YOLD,Y) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION Y(N),YOLD(N),DUMMY(N),FUNOLD1(N),FUNOLD2(N) COMMON x1,x2,h,phi1,phi2 c istep = 0 Time = x1 c c Print the Initial solution to output file write(11,*) write(11,'(12x,a,13x,a,13x,a)') 'X','Y(1)','Y(2)' write(11,'(3f16.6)') Time,YOLD(1),YOLD(2) c c Start Integration Loop 10 Time = Time + h istep = istep + 1 c c Set Termination Criteria for the Integration Loop if((Time-x2) .gt. 1.d-14) return c c Evaluate the RHS Function at Previous Time call Func(N,YOLD,FUNOLD1) c c For the first time step apply 1st-order explicit Euler for predictor c and 1st-order implicit Euler for corrector if (istep .eq. 1) then c predicor do i = 1, N Y(i) = YOLD(i)+h*FUNOLD1(i) enddo c evaluate RHS function at predicted variable call Func(N,Y,DUMMY) c corrector do i = 1, N Y(i) = YOLD(i)+h*DUMMY(i) enddo else c c Apply 2nd-order predictor-corrector method for Remaining Time Steps c predicor do i = 1, N Y(i) = YOLD(i)+h/2.d0*(3.d0*FUNOLD1(i)-FUNOLD2(i)) enddo c evaluate RHS function at predicted variable call Func(N,Y,DUMMY) c corrector do i = 1, N Y(i) = YOLD(i)+h/2.d0*(FUNOLD1(i)+DUMMY(i)) enddo endif c c Print the solution for the current time step to output file write(11,'(3f16.6)') Time,Y(1),Y(2) c c Store Previous Time Values do i = 1, N YOLD(i) = Y(i) FUNOLD2(i) = FUNOLD1(i) enddo go to 10 c end c c c c Subroutine Adams_Moulten to integrate IVP's c c c subroutine Adams_Moulten(N,ITMAX,TOL,IFLAG,YOLD,Y) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION Y(N),Y_G(N),FUN(N),FUNOLD(N),YOLD(N) COMMON x1,x2,h,phi1,phi2 c istep = 0 Time = x1 c c Print the Initial solution to output file write(11,*) write(11,'(12x,a,13x,a,13x,a)') 'X','Y(1)','Y(2)' write(11,'(3f16.6)') Time,YOLD(1),YOLD(2) c c Start Integration Loop 10 Time = Time + h istep = istep + 1 write(*,'(4x,a,f10.4)') 'Time =',Time c c Set Termination Criteria for the Integration Loop if((Time-x2) .gt. 1.d-14) return c c Evaluate the RHS Function for previous time call Func(N,YOLD,FUNOLD) c c Initial guess for current integration step is the solution form previous step do i = 1, N Y_G(i) = YOLD(i) enddo call Newton(N,ITMAX,TOL,IFLAG,ISTEP,YOLD,FUNOLD,Y_G,Y) c c Evaluate the RHS Function for current time call Func(N,Y,FUN) c c Store previous integration values do i = 1, N YOLD(i) = Y(i) FUNOLD(i) = FUN(i) enddo c c Print the solution for the current time step to output file write(11,'(3f16.6)') Time,Y(1),Y(2) c go to 10 c return end c c c c Subroutine RKG4 to integrate IVP's using Rung Kutta Gill 4th order method c c c subroutine RKG4(N,YOLD,Y) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION Y(N),FUN(N),FUNOLD(N),YOLD(N) DOUBLE PRECISION rk1(N),rk2(N),rk3(N),rk4(N),Y_DUMMY(N),F_DUMMY(N) COMMON x1,x2,h,phi1,phi2 c a = (dsqrt(2.d0)-1.d0)/2.d0 b = (2.d0-dsqrt(2.d0))/2.d0 c = -dsqrt(2.d0)/2.d0 d = 1.d0+dsqrt(2.d0)/2.d0 c istep = 0 Time = x1 c c Print the Initial solution to ivp.txt write(11,*) write(11,'(12x,a,13x,a,13x,a)') 'X','Y(1)','Y(2)' write(11,'(3f16.6)') Time,YOLD(1),YOLD(2) c c Evaluate the RHS Function for the Initial Solution call Func(N,YOLD,FUNOLD) c c Start Integration Loop 10 Time = Time + h istep = istep + 1 c c Set Termination Criteria for the Integration Loop if((Time-x2) .gt. 1.d-14) return c do i = 1, N rk1(i) = h*FUNOLD(i) Y_DUMMY(i) = YOLD(i)+0.5d0*rk1(i) enddo call Func(N,Y_DUMMY,F_DUMMY) c do i = 1, N rk2(i) = h*F_DUMMY(i) Y_DUMMY(i) = YOLD(i)+a*rk1(i)+b*rk2(i) enddo call Func(N,Y_DUMMY,F_DUMMY) c do i = 1, N rk3(i) = h*F_DUMMY(i) Y_DUMMY(i) = YOLD(i)+c*rk2(i)+d*rk3(i) enddo call Func(N,Y_DUMMY,F_DUMMY) c do i = 1, N rk4(i) = h*F_DUMMY(i) enddo c do i = 1, N Y(i) = YOLD(i)+1.d0/6.d0*(rk1(i)+rk4(i))+ ! 1.d0/3.d0*(b*rk2(i)+d*rk3(i)) enddo c call Func(N,Y,FUN) c c Store previous integration values do i = 1, N YOLD(i) = Y(i) FUNOLD(i) = FUN(i) enddo c c Print the solution for the current time step to ivp.txt write(11,'(3f16.6)') Time,Y(1),Y(2) c go to 10 c return end c c c c c Subroutine Newton to to solve a system of algebraic equations iteratively c c c subroutine Newton !(N,ITMAX,TOL,IFLAG,ISTEP,YOLD,FUNOLD,UNKNON_G,UNKNON) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION A(N,N),UNKNON_G(N),UNKNON(N),B(N) DOUBLE PRECISION YOLD(N),FUNOLD(N) integer ipvt(N) COMMON x1,x2,h,phi1,phi2 c c initialize matrix A and vector B do i = 1, N UNKNON(i) = UNKNON_G(i) 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(N,YOLD,FUNOLD,UNKNON,B) c c Check the norm of the residuals rnorm = 0.d0 do i = 1, N rnorm = rnorm + B(i)**2 enddo rnorm = dsqrt(rnorm) write(*,'(3x,a,i3,5x,a,e20.10)') ! 'iteration #',iter, 'Tolerance =', rnorm if (rnorm .lt. tol) return c if (istep .gt. 10) go to 20 c c Evaluate the Jacobian Matrix call jacobian(IFLAG,N,YOLD,FUNOLD,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,N,N,ipvt,d) 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 20 do i = 1, N B(i) = -B(i) enddo c call lubksb(A,N,N,ipvt,B) c do i= 1, N 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 Subroutine Residuals to Input System of Non-Linear Equations c c c subroutine Residuals(N,YOLD,FUNOLD,UNKNON,R) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION UNKNON(N),YOLD(N),FUN(N),FUNOLD(N),R(N) COMMON x1,x2,h,phi1,phi2 c do i = 1, N R(i) = 0.d0 enddo c c Evaluate the RHS Functions of IVP's for UNKNON call Func(N,UNKNON,FUN) c c Evaluate the Residual Equations do i = 1, N R(i) = UNKNON(i)-YOLD(i)-h/2.d0*(FUN(i)+FUNOLD(i)) enddo c return end c c c c Subroutine jacobian to Evaluate the Jacobian Matrix c c c subroutine jacobian(IFLAG,N,YOLD,FUNOLD,UNKNON,A) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION A(N,N),UNKNON(N),UNKNONJ(N),FUN(N),FUNJ(N) DOUBLE PRECISION YOLD(N),FUNOLD(N) COMMON x1,x2,h,phi1,phi2 parameter (epsl=1.d-5) c do i = 1, N do j = 1, N A(i,j)=0.d0 enddo enddo c if (IFLAG .gt. 0) go to 10 c c Evaluate the jacobian matrix numerically 10 do i = 1, N UNKNONJ(i)=UNKNON(i) enddo call residuals(N,YOLD,FUNOLD,UNKNON,FUN) do i = 1, N diff = dmax1(epsl,dabs(epsl*UNKNON(i))) UNKNONJ(i) = UNKNON(i)+diff call residuals(N,YOLD,FUNOLD,UNKNONJ,FUNJ) do j = 1, N A(j,i) = (FUNJ(j)-FUN(j))/diff enddo UNKNONJ(i)=UNKNON(i) enddo 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