C*********************************************************************** C * C BROYDEN'S METHOD * C * C*********************************************************************** C C C C TO APPROXIMATE THE SOLUTION OF THE NONLINEAR SYSTEM F(X)=0 C GIVEN AN INITIAL APPROXIMATION X. C C INPUT: NUMBER N OF EQUATIONS AND UNKNOWNS; INITIAL C APPROXIMATION X=(X(1),...,X(N)); TOLERANCE TOL; C MAXIMUM NUMBER OF ITERATIONS N. C C OUTPUT: APPROXIMATE SOLUTION X=(X(1),...,X(N)) OR A MESSAGE C THAT THE NUMBER OF ITERATIONS WAS EXCEEDED. C DIMENSION A(3,3),C(3,3),X(3),Y(3),S(3),Z(3),V(3) C USE NN INSTEAD OF N FOR THE MAXIMUM NUMBER OF ITERATIONS CHARACTER NAME*14,NAME1*14,AAA*1 INTEGER INP,OUP,FLAG LOGICAL OK OPEN(UNIT=5,FILE='CON',ACCESS='SEQUENTIAL') OPEN(UNIT=6,FILE='CON',ACCESS='SEQUENTIAL') WRITE(6,*) 'This is the Broyden Method for Nonlinear Systems.' WRITE(6,*) 'Has the function F been defined and have the' WRITE(6,*) 'partial derivatives been defined as follows:' WRITE(6,*) '1. function F(N,I,X) where I is the number of the' WRITE(6,*) ' component function, X is the vector at which' WRITE(6,*) ' to evaluate and N is the dimension.' WRITE(6,*) '2. function FP(N,I,J,X) where I is the number of' WRITE(6,*) ' the component function and J is the number' WRITE(6,*) ' of the variable with respect to which' WRITE(6,*) ' partial differentiation is performed,' WRITE(6,*) ' X is the vector at which to evaluate and N' WRITE(6,*) ' is the dimension.' WRITE(6,*) 'Enter Y or N ' WRITE(6,*) ' ' READ(5,*) AAA IF(( AAA .EQ. 'Y' ) .OR. ( AAA .EQ. 'y' )) THEN OK = .FALSE. 101 IF (OK) GOTO 11 WRITE(6,*) 'Input the number n of equations.' WRITE(6,*) ' ' READ(5,*) N IF (N.GE.2) THEN OK = .TRUE. ELSE WRITE(6,*) 'N must be greater than 1.' ENDIF GOTO 101 11 OK = .FALSE. 14 IF (OK) GOTO 15 WRITE(6,*) 'Input tolerance ' WRITE(6,*) ' ' READ(5,*) TOL IF (TOL.LE.0.0) THEN WRITE(6,*) 'Tolerance must be positive ' ELSE OK = .TRUE. ENDIF GOTO 14 15 OK=.FALSE. 12 IF (OK) GOTO 13 WRITE(6,*) 'Input the maximum number of iterations.' WRITE(6,*) ' ' READ(5,*) NN IF ( NN .LE. 0 ) THEN WRITE(6,*) 'Must be positive integer ' ELSE OK = .TRUE. ENDIF GOTO 12 13 DO 16 I=1,N WRITE(6,*) 'Input initial approximation X(',I,')' WRITE(6,*) ' ' READ(5,*) X(I) 16 CONTINUE ELSE WRITE(6,*) 'The program will end so that the functions' WRITE(6,*) 'F(N,I,X) and FP(N,I,J,X) can be created ' OK = .FALSE. ENDIF IF(.NOT.OK) GOTO 400 WRITE(6,*) 'Select output destinations: ' WRITE(6,*) '1. Screen ' WRITE(6,*) '2. Text file ' WRITE(6,*) 'Enter 1 or 2 ' WRITE(6,*) ' ' READ(5,*) FLAG IF ( FLAG .EQ. 2 ) THEN WRITE(6,*) 'Input the file name in the form - ' WRITE(6,*) 'drive:name.ext' WRITE(6,*) 'with the name contained within quotes' WRITE(6,*) 'as example: ''A:OUTPUT.DTA'' ' WRITE(6,*) ' ' READ(5,*) NAME1 OUP = 3 OPEN(UNIT=OUP,FILE=NAME1,STATUS='NEW') ELSE OUP = 6 ENDIF WRITE(OUP,*) 'BROYDEN METHOD FOR NONLINEAR SYSTEMS' WRITE(OUP,*) 'Iteration, approximation x(1),...,x(n) and error' PI = 4*ATAN(1.) C SN WILL REPRESENT THE L2 NORM OF THE ERROR SN=0 K=0 WRITE(OUP,3) K,(X(I),I=1,N),SN C STEP 1 C A WILL HOLD THE JACOBIAN FOR THE INITIAL APPROXIMATION C THE SUBPROGRAM FP COMPUTES THE ENTRIES OF THE JACOBIAN DO 10 I=1,N DO 10 J=1,N 10 A(I,J)=FP(N,I,J,X) C COMPUTE V = F(X(0)) C THE SUBPROGRAM F(N,I,X) COMPUTES THE ITH COMPONENT OF F AT X DO 20 I=1,N 20 V(I)=F(N,I,X) C STEP 2 CALL INVER(N,A,C) C INVER FINDS THE INVERSE OF THE N BY N MATRIX A AND RETURNS IT IN A C C IS A DUMMY PARAMETER USED TO RESERVE STORAGE IN THE SUBROUTINE C STEP 3 K=1 C NOTE: S=S(1) CALL MULT(N,A,V,S,SN) C MULT COMPUTES THE PRODUCT S=-A*V AND THE L2-NORM SN OF S DO 30 I=1,N 30 X(I)=X(I)+S(I) C NOTE: X=X(1) WRITE(OUP,3) K,(X(I),I=1,N),SN C STEP 4 200 IF ( K.GT.NN ) GOTO 300 C STEP 5 C THE VECTOR W IS NOT USED SINCE THE FUNCTION F IS EVALUATED C COMPONENT BY COMPONENT DO 40 I=1,N VV=F(N,I,X) Y(I)=VV-V(I) 40 V(I)=VV C NOTE: V=F(X(K)) AND Y=Y(K) C STEP 6 CALL MULT(N,A,Y,Z,ZN) C NOTE: Z=-A(K-1)**-1 * Y(K) C STEP 7 P=0 C P WILL BE S(K)**T * A(K-1)**-1 * Y(K) DO 60 I=1,N P=P-S(I)*Z(I) C STEP 8 DO 60 J=1,N 60 C(I,J)=(S(I)+Z(I))*S(J) DO 100 I=1,N 100 C(I,I)=C(I,I)+P C C=S(K)**T * A(K-1)**-1 * Y(K)*I + (S(K)+A(K-1)**-1 * Y(K))* C S(K)**T ) C STEPS 9 AND 10 DO 70 I=1,N DO 80 J=1,N ACC=0 DO 90 L=1,N C ACC IS THE (I,J) ENTRY OF THE INVERSE OF A 90 ACC=ACC+C(I,L)*A(L,J)/P C Z WILL HOLD THE J-TH COLUMN OF THE INVERSE OF A 80 Z(J)=ACC DO 70 J=1,N 70 A(I,J)=Z(J) CALL MULT(N,A,V,S,SN) C NOTE: A=A(K)**-1 AND S=-A(K)**-1 * F(X(K)) C STEP 11 DO 110 I=1,N 110 X(I)=X(I)+S(I) C NOTE: X = X(K+1) WRITE(OUP,3) K+1,(X(I),I=1,N),SN C STEP 12 IF (SN.LE.TOL) THEN C PROCEDURE COMPLETED SUCCESSFULLY WRITE(OUP,6) GOTO 400 END IF C STEP 13 K=K+1 GOTO 200 C STEP 14 300 WRITE(OUP,4) 400 CLOSE(UNIT=5) CLOSE(UNIT=OUP) IF(OUP.NE.6) CLOSE(UNIT=6) STOP 3 FORMAT(1X,I2,3(1X,E15.8,1X),E15.8) 4 FORMAT(1X,'MAXIMUM NUMBER OF ITERATIONS EXCEEDED') 6 FORMAT(1X,'PROCEDURE COMPLETED SUCCESSFULLY') END C FUNCTION FP(N,I,J,X) DIMENSION X(N) IF (I.EQ.1) THEN IF(J.EQ.1) FP=3 IF(J.EQ.2) FP=X(3)*SIN(X(2)*X(3)) IF(J.EQ.3) FP=X(2)*SIN(X(2)*X(3)) ELSE IF (I.EQ.2) THEN IF(J.EQ.1) FP=2*X(1) IF(J.EQ.2) FP=-162*(X(2)+.1) IF(J.EQ.3) FP= COS(X(3)) ELSE IF(J.EQ.1) FP=-X(2)*EXP(-X(1)*X(2)) IF(J.EQ.2) FP= -X(1)*EXP(-X(1)*X(2)) IF(J.EQ.3) FP=20 ENDIF ENDIF RETURN END C FUNCTION F(N,I,X) DIMENSION X(N) IF(I.EQ.1) F=3*X(1)-COS(X(2)*X(3))-.5 IF(I.EQ.2) F=X(1)**2-81*(X(2)+.1)**2+SIN(X(3))+1.06 IF(I.EQ.3) F=EXP(-X(1)*X(2))+20*X(3)+(10*3.141593-3)/3 RETURN END C SUBROUTINE MULT(N,A,Y,Z,ZN) DIMENSION A(N,N),Y(N),Z(N) ZN=0 DO 10 I=1,N Z(I)=0 DO 15 J=1,N 15 Z(I)=Z(I)-A(I,J)*Y(J) 10 ZN=ZN+Z(I)*Z(I) ZN=SQRT(ZN) RETURN END C SUBROUTINE INVER(N,A,B) DIMENSION A(N,N),B(N,N) DO 10 I=1,N DO 10 J=1,N B(I,J)=0 10 IF(J.EQ.I) B(I,J)=1 DO 20 I=1,N I1=I+1 I2=I IF(I.NE.N) THEN DO 30 J=I1,N 30 IF(ABS(A(J,I)).GT.ABS(A(I2,I))) I2=J IF(I2.NE.I) THEN DO 40 J=1,N C=A(I,J) A(I,J)=A(I2,J) A(I2,J)=C C=B(I,J) B(I,J)=B(I2,J) B(I2,J)=C 40 CONTINUE END IF END IF DO 50 J=1,N IF (J.NE.I) THEN C=A(J,I)/A(I,I) DO 60 K=1,N A(J,K)=A(J,K)-C*A(I,K) B(J,K)=B(J,K)-C*B(I,K) 60 CONTINUE END IF 50 CONTINUE 20 CONTINUE DO 70 I=1,N C=A(I,I) DO 70 J=1,N 70 A(I,J)=B(I,J)/C RETURN END