c c Sample_Code_2 c Solves a System of Non-Linear Equations by Newton-Raphson Method c R(X) = 0 c IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (N = 2) DOUBLE PRECISION UNKNON_G(N),UNKNON(N) OPEN(11,FILE='sample_code_2_output.txt') c itmax = 10 ! Maximum Number of Iterations tol = 1.d-10 ! Maximum Tolerance iflag = 1 ! Flag for the Jacobian Matrix (0 -> analytical Jacobian, 1 -> Numerical Jacobian) c c Provide an initial guess UNKNON UNKNON_G(1) = 1.5d0 UNKNON_G(2) = 3.5d0 c call Newton to solve the non-linear equations iteratively call Newton(N,itmax,tol,iflag,UNKNON_G,UNKNON) c c write the solution to newton.txt 11 do i = 1, N write(* ,'(e20.10)') UNKNON(i) write(11,'(e20.10)') UNKNON(i) enddo end c c c c c c c c Subroutine Newton to to solve a system of algebraic equations iteratively subroutine Newton(N,itmax,tol,iflag,UNKNON_G,UNKNON) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION A(N,N),UNKNON_G(N),UNKNON(N),B(N) integer ipvt(N) c c initialize matrix A and vector B do i = 1, N UNKNON(i) = UNKNON_G(i) B(i) = 0.0 do j = 1, N 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(N,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 c Evaluate the Jacobian Matrix c IFLAG = 0 ==> analatical jacobian c IFLAG = 1 ==> numerical jacobian call jacobian(IFLAG,N,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 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 c c c c c c Subroutine Residuals to Input System of Non-Linear Equations subroutine residuals(N,UNKNON,R) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION UNKNON(N),R(N) c do i = 1, N R(i) = 0.d0 enddo c R(1) = UNKNON(1)**2+UNKNON(1)*UNKNON(2)-10.d0 R(2) = UNKNON(2)+3.d0*UNKNON(1)*UNKNON(2)**2-57.d0 c return end c c c c c c c c c Subroutine jacobian to Evaluate the Jacobian Matrix subroutine jacobian(IFLAG,N,UNKNON,A) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION A(N,N),UNKNON(N),UNKNONJ(N),FUN(N),FUNJ(N) 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 analatically A(1,1) = 2.d0*UNKNON(1)+UNKNON(2) A(1,2) = UNKNON(1) A(2,1) = 3.d0*UNKNON(2)**2 A(2,2) = 1.d0+6.d0*UNKNON(1)*UNKNON(2) return c c Evaluate the jacobian matrix numerically 10 do i = 1, N UNKNONJ(i)=UNKNON(i) enddo call residuals(N,UNKNON,FUN) do i = 1, N diff = dmax1(epsl,dabs(epsl*UNKNON(i))) UNKNONJ(i) = UNKNON(i)+diff call residuals(N,UNKNONJ,FUNJ) do j = 1, N A(j,i) = (FUNJ(j)-FUN(j))/diff enddo UNKNONJ(i)=UNKNON(i) enddo return end 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