LIBRARY OF THE 
 
 UNIVERSITY OF ILLINOIS 
 
 AT URBANA-CHAMPAIGN 
 
 5IO.ft4 
 
The person charging this material is re- 
 sponsible for its return to the library from 
 which it was withdrawn on or before the 
 Latest Date stamped below. 
 
 Theft, mutilation, and underlining of books 
 are reasons for disciplinary action and may 
 result in dismissal from the University. 
 
 UNIVERSITY OF ILLINOIS LIBRARY AT URBANA-CHAMPAIGN 
 
 DEC 1 4 1976 
 
 DEC 12^ 
 
 «*10 
 
 y. 
 
 xii 
 
 AUG 2 9 1997 
 
 L161 — O-1096 
 
J/"' 
 
 UIUCDCS-R-73-576 
 
 77? (Ltt, 
 
 coo-ll+69-0226 
 
 
 EVALUATION OF AN INTEGRATION TECHNIQUE FOR THE 
 SOLUTION OF NONLINEAR EQUATIONS 
 
 "by 
 
 James Harold Robb 
 
 August 19 73 
 
 DEPARTMENT OF COMPUTER SCIENCE 
 UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 
 
 URBANA, ILLINOIS 
 
 THE LIBRARY OF THE 
 
 ^'y.ERSITY OF ILLINOIS 
 AT URRAMa <"""mdaigN 
 
UIUCDCS-B-73-576 
 
 EVALUATION OF AN INTEGRATION TECHNIQUE FOR THE 
 SOLUTION OF NONLINEAR EQUATIONS* 
 
 by 
 James Harold Robb 
 
 August 1973 
 
 DEPARTMENT OF COMPUTER SCIENCE 
 
 UNIVERSITY OF ILLINOIS AT URBAN A- CHAMPAIGN 
 
 URBANA, ILLINOIS 6l801 
 
 * Supported in part by the Atomic Energy Commission under contract 
 US AEC AT(ll-l)lU69 and submitted in partial fulfillment of the 
 requirements of the Graduate College for the degree of Master of 
 Science in Computer Science. 
 
Digitized by the Internet Archive 
 in 2013 
 
 http://archive.org/details/evaluationofinte576robb 
 
Ill 
 
 ACKNOWLEDGMENT 
 
 The author is grateful for the guidance and suggestions 
 received in a number of helpful discussions with his advisor, 
 Professor C. W. Gear, and expresses his thanks for that help. Thanks 
 are also given to Barbara Armstrong for her very efficient work in 
 typing the thesis. 
 
IV 
 
 TABLE OF CONTENTS 
 
 Page 
 
 1. INTRODUCTION 1 
 
 2. SUBROUTINE INTECH 3 
 
 2.1 Description of the Method 3 
 
 2.2 Strategy 5 
 
 2.3 The Calling Program and Other Subroutines 8 
 
 3. NUMERICAL RESULTS 10 
 
 3.1 Description of the Examples 11 
 
 3.2 Discussion of Results 15 
 
 k. CONCLUSIONS 22 
 
 LIST OF REFERENCES 2.k 
 
 APPENDIX 26 
 
1. INTRODUCTION 
 
 The method for solving systems of nonlinear algebraic equations 
 that is the subject of this investigation was developed by van Me lie [l]. 
 His subroutine DIFMF3 was designed to solve the steady-state problem of 
 ordinary differential equations , and for this purpose each component of 
 the vector to be solved for could be either a derivative or an ordinary 
 variable. The development of this method grew out of ideas presented by 
 Gear [2,3] concerning the solution of systems of mixed algebraic and 
 differential equations. 
 
 In order to evaluate this method simply as a means of solving 
 nonlinear equations , another subroutine called INTECH was prepared by 
 the author. This subroutine is essentially the same as DIFMF3 but is 
 simplified so as to handle only a single vector of variables instead of 
 giving the option of choosing each component from either a vector of 
 variables or a vector of derivatives. Provisions for applying constraints 
 to the variables were also omitted, along with certain variable names that 
 were used in DIFMF3 only to pass data to other subroutines in the 
 differential -equation-solving package of which it was a part. Finally, 
 certain changes were required to convert the subroutine from FORTRAN IV 
 to CALL OS FORTRAN, which the author chose to use because of the 
 convenience of the timesharing system available for this language. 
 
 Most of the examples used in the investigation were taken from 
 articles by Boggs [U,5] or by Powell [6,7] so that comparisons could be 
 made with their published results. Their results appeared to the author 
 to be good ones to compare with because they both report on recent 
 developments, their methods are different, and both show some improvement 
 
. _„__! ■'_.-"* 
 
 
2. SUBROUTINE INTECH 
 
 This section presents a description of the "basic method employed, 
 the strategy used in an attempt to economize on the amount of computation 
 required to obtain a solution, and a discussion of some of the details to 
 be considered in using the subroutine. For a listing of the subroutine, 
 along with the main program and subroutines used for one of the examples, 
 refer to the Appendix. 
 
 2. 1 Description of the Method 
 
 The problem that we wish to solve may be written in the form 
 
 f(y) = (2.1) 
 
 r 1 2 n i T 
 
 where y = Ly , y >• • ■ y J 
 
 12 _n T 
 
 and f(y) = [f (y), f (y),... i (y)] are n-component column vectors, and 
 
 is an n-component column vector of zeros. Also, we will use J(y) to 
 
 denote the Jacobian matrix 
 
 J(y) -&&k . 
 
 y 
 
 Following Boggs ' approach in Section 2 of his paper [U], we consider the 
 operator homotopy G(s,y) = which imbeds the real parameter s into the 
 original equation so that G(0,y) = has the solution y , and 
 G(s,y) ► f(y) as s ■*- °°. Specifically, we use 
 
 G(s,y) = f(y) - e" S f(y ) = 
 to achieve these results. Then, since G(s,y) = 0, we may obtain y' 
 
 dy_ 
 ds 
 
 rale , in the form 
 
 (denoting r~- ) by differentiating with respect to s and using the chain 
 
 y' = -J 1 (y)f(y), y(0) = y Q . (2.2) 
 
Euler's method of numerical integration applied to this differential 
 equation gives the iteration 
 
 y m = y m -i + hy ;-i (2 - 3) 
 
 We note that this is equivalent, if we make h = 1, to Newton's method 
 applied to the original equation (2.1). Formula (2.3) is used as the 
 predictor step in van Melle's method. For the corrector, he uses 
 
 y m = y m-l +h(ay m + (l " a) K-l ] {2 - k) 
 
 with <_ a £ 1. When a = 0, (2.U) is the same as (2.3). If a = 
 
 and h = 1, the iteration of the predictor and corrector steps is identical 
 
 to Newton's method applied to (2.1), and if it converges to a solution, 
 
 will do so relatively quickly. On the other hand, if a = 1 and h << 1, 
 
 the method is no longer Newton's method and will generally require many 
 
 steps to converge, but the convergence is more likely to occur. 
 
 Note that, as it stands, (2.1+) cannot be solved directly for 
 
 y because y', which is a function of y by (2.2), appears on the right- 
 m m m 
 
 hand side. Rather than solving (2.k) iteratively each time, thus getting 
 
 an implicit method, we obtain an approximate solution in an explicit 
 
 form. To do this, denote the right-hand side of (2.3), which is the 
 
 result of the predictor operation, by y , then denote the right-hand 
 
 side of (2.1+) by y = y + Ay . Using a first order Taylor expansion 
 
 mm m 
 
 on the right-hand side of (2.2) to represent y' , the following explicit 
 corrector may be obtained so that the method becomes 
 
A ■ est t hJ_1 < 0> > + ^-J < 2 - 6 » 
 
 (O) , V 
 
 y m =y m - «a (2.7) 
 
 hy ; = hy ;.i- A (2 ' 8) 
 
 where (2.5) is the predictor and (2.6) through (2.8) make up the 
 
 corrector. A single iteration of the corrector is used after each 
 
 application of the predictor. 
 
 If the inverse Jacobian were to "be evaluated at every step, 
 
 we would use J (y ) in (2.6). However, part of the strategy in 
 m 
 
 applying this method is an attempt to reduce the amount of computation 
 required by evaluating the Jacobian and its inverse infrequently, as 
 discussed in Section 2.2. Therefore, the symbol J is used in (2.6) 
 to denote the inverse Jacobian that is currently being used. 
 
 2.2 Strategy 
 
 One of the major problems that had to be solved empirically 
 was the method of varying a and h in order to find a solution as 
 efficiently as possible. The method used in INTECH is basically the 
 same as that used in van Me lie's subroutine DIFMF3 [l], with any 
 differences being noted in the description that follows. In order to 
 relate the discussion more closely to the listing of the subroutine 
 INTECH in the Appendix, the variable names ALPHA and H will be used to 
 denote, respectively, the symbols a and h used in Section 2.1. 
 
Initially, the procedure begins with ALPHA = and H = 1, 
 so that the method is the same as Newton's method. As long as the norm 
 
 l|f(y)M = * ly 1 ! 
 
 i=l 
 decreases significantly, ALPHA and H are not changed. The norm is 
 calculated following each predictor step. If the norm fails to decrease 
 by at least five percent from its previous value, the replacements 
 ALPHA = 1 and H = .01 are made and the previous vector y is restored 
 before repeating the step. After this, on each succeeding iteration as 
 long as the norm has decreased by at least two percent from its previous 
 value, ALPHA is decreased by multiplying it by the factor 0.8, and H is 
 increased by multiplying it by the factor R = 1.7 - . 85H + . 15/H which 
 is designed so that the limitation H <_ 1 will be observed. This is one 
 point at which INTECH differs from van Me lie's subroutine DIFMF3. In 
 DIFMF3 the corresponding formula for R is: R = .6 + .k/E. In preliminary 
 tests on INTECH involving just Examples 1 through 9» it was found that 
 when using the latter formula for R, no convergence was obtained for 
 Examples 5 and 9- The indicated formula for R was found to be the best 
 of several that were tried, giving convergence for all nine examples and 
 requiring the same or a smaller number of function evaluations than the 
 original formula for all except Example 3 9 in which case one additional 
 evaluation was required. 
 
 If on any iteration after the Initially used Newton's method 
 has failed, the norm exceeds the previous norm by a factor of 100 or more, 
 the replacements ALPHA = 1 and H = max(H/2, 0.2) are made, the previous 
 vector y is restored, and the predictor step is repeated. However, if the 
 
7 
 
 ratio of the norm to the previous norm remains between O.98 and 100, 
 the replacements ALPHA = 1 and R = min(l.3, . 6/H) are made and the 
 corrector step is performed. 
 
 Another feature of the program related to strategy is the 
 separation of the vector y of variables into two vectors, given the 
 names Y and YL in the subroutine INTECH. The values of all variables 
 that appear nonlinear ly in one or more of the functions are stored in 
 Y; while those that appear only linearly are stored in YL. This 
 separation makes it convenient to take advantage of the fact, which 
 is discussed by Gear [2,3], that only the correction process and not 
 the prediction process needs to be applied to the linear variables. 
 For systems involving a large number of linear variables , this feature 
 may reduce the amount of computation considerably. 
 
 Another problem that had to be considered was the frequency of 
 evaluating and inverting the Jacobian. It is not desirable to do this 
 on every step, particularly where n is large, because of the excessive 
 amount of computation required. The strategy that has been developed 
 calls for recomputing the Jacobian every 5NY steps where NY is the 
 number of variables in the vector Y discussed in the preceding paragraph. 
 However, if the Jacobian has not been evaluated recently and, in 
 addition, either the norm increases by a factor of 100 or more from its 
 previous value or else it becomes less than 1.0 (indicating that a 
 solution is being approached) , then the Jacobian is computed and inverted 
 even though this action is not called for by the first rule. 
 
8 
 
 2 .3 The Calling Program and Other Subroutines 
 
 The calling program must define the variables required by the 
 parameter list, which is given in the listing of INTECH in the Appendix. 
 Appropriate input values must be supplied for NSEND, the maximum number 
 of iterations to be performed; DEL, the value used to compare with the 
 norm 
 
 n 
 
 * |y x l 
 
 i=l 
 to determine whether or not satisfactory convergence has been obtained; 
 NY, the number of nonlinear variables (the dimension of Y); NL , the 
 number of linear variables (the dimension of YL); N, the total number 
 of variables, NY + NL ; NY starting values in the array Y, and NL 
 starting values in the array YL. Results are printed by statements 
 within INTECH after each prediction process of the iteration, so no 
 provisions for output are normally needed in the calling program. The 
 formats used in INTECH were designed for use with a terminal taking 
 paper 8 1/2" wide, and should be changed appropriately if a line 
 printer is to be used. 
 
 In addition to the calling program, the user must supply 
 two subroutines that are tailored to the specific problem, called 
 DIFFUN and MATSET, for which examples are given in the Appendix. The 
 purpose of DIFFUN is to calculate the vector f with the value of f 
 being stored in DY(l). The Jacobian matrix is calculated by MATSET 
 
 and stored in the array PW. On return from MATSET, the contents of 
 
 <-) tw ( T \ 
 PW(I,J) are ( — T ) for I from 1 to N and J from 1 to NY, and 
 
 a I { J ) 
 
 , 3DY ( I ) N 
 9YL(J NY) ' f ° r fr ° m 1 t0 N and J from WY+1 to N * 
 
Finally, two more subroutines are required that do not 
 depend on the specific problem to be solved. The first of these 
 is MINV, a matrix inversion subroutine that takes as input the 
 Jacobian matrix in the array PW and provides as output the inverse 
 Jacobian in the same array. Since a number of matrix inversion 
 routines that are generally available may be adapted, the particular 
 subroutine used in this investigation is not given in the Appendix. 
 The second of these two is MATMUL, a short subroutine that multiplies 
 the function vector in DY by the inverse Jacobian in PW with the 
 resulting vector being stored in Fl on return. A listing of this 
 subroutine is included in the Appendix. 
 
10 
 3. NUMERICAL RESULTS 
 
 In order to have a good basis for comparing the numerical results 
 from INTECH with the results obtained from programs using other methods for 
 the solution of nonlinear equations, it was decided to choose most of the 
 examples from previously published and widely available investigations made 
 by others, and to make comparisons with these published results. This 
 avoids the possibility of imperfectly implementing the various authors' 
 ideas into computer programs, which could lead to invalid conclusions. 
 However, it does not make it possible to have the comparisons as straight- 
 forward as would be desirable, in some cases, because of different 
 approaches and different kinds and amounts of detail being presented by 
 different authors. 
 
 In addition to making comparisons with the published results 
 of others, it was decided to also run each problem on an already available 
 library subroutine for finding the minimum of a function based on the 
 Fletcher-Powell algorithm [8] as detailed by Wells [9] and later modified 
 by Fletcher [10]. This involved only making modifications to translate 
 the available subroutine from FORTRAN IV to CALL OS FORTRAN so that it 
 could be used in the timesharing mode along with INTECH. This modified 
 Fletcher-Powell subroutine is referred to as FLPMIN in the discussion of 
 results. The function to be minimized by this subroutine was taken to 
 be the sum of the squares of the functions that formed the system to be 
 solved simultaneously in each case. 
 
11 
 
 3.1 Description of the Examples 
 
 Example 1 
 
 Function: 
 
 f(x) = 
 
 h + x, + x_ - x. + 2x n x + 3x, 
 
 1 2 
 
 1 + 2x ± - 3x 2 + x ± + x 1 x 2 - 2x, 
 
 Starting Vector: x = (-2.057, -7-503) 
 Approximate Solution: x= (3-339, -2.98*0 
 
 This example is one of those used by van Melle [l] in evaluating 
 the subroutine DIFMF3 upon which INTECH is based. 
 
 Example 2 
 
 Function: same as Example 1 
 
 Starting Vector: x = (0,1) 
 
 Approximate Solution: x = (-1.533U, 0.06ll2l) 
 
 It was noted that Example 1 had another solution, and this 
 starting vector was tried with the expectation of obtaining the 
 above solution. 
 
 Example 3 
 Function: 
 
 f(x) = 
 
 2 
 
 X l " 
 
 *2 
 
 + 1 
 
 x n - cos (J- x_) 
 
 1 _x 2 *2 
 
 Starting Vector: x = (1,0) 
 Solution: x = (0,1) 
 
 This example is discussed in detail by Boggs [U,5]. The 
 solution given is the one which can be located from this starting point 
 without crossing the lines of singularity for the Jacobian. 
 
12 
 
 Example k 
 
 Function: same as Example 3 
 Starting Vector: x = (-1,1) 
 Solution: (- — , 1.5) 
 
 This particular starting vector was not discussed by Boggs, 
 but it is representative of the starting region which he states should 
 lead to the given solution [H,5]. 
 
 Example 5 
 Function: 
 
 f(x) = 
 
 1 *2 
 
 - [sinU^) - — - x 1 ] 
 
 t 1 " 57> t 
 
 2x ± ey^ 
 
 e - e] + — - 2ex 1 
 
 Starting Vector: x = (O.U, 3) 
 Approximate Solution: x = (0.3, 2.8) 
 
 This is the second of the two examples discussed by Boggs [4,5] 
 and it is also given in an earlier article by Broyden [ll]. 
 
 Example 6 
 Fun ct i on : 
 
 f(x) = 
 
 10x l 2 
 
 Z^T +2x 2 
 
 Starting Vector: x = (3,1) 
 Solution: x = (0,0) 
 
 This is one of the examples discussed by Powell [6], One of his 
 earlier methods did not converge to the correct solution of this problem, 
 owing to the singularity of the Jacobian matrix at points between the starting 
 
13 
 
 point and the solution. Powell indicated that his new algorithm was 
 developed in order to provide convergence for examples such as this, 
 but information on the number of function evaluations required for 
 this particular example is apparently not given in his discussion of 
 numerical results [7]. 
 
 Example 7 
 Function: 
 
 f(x) = 
 
 IOjOOOx^x^ 
 
 - 1 
 
 -x^ 
 
 -X 1 "2 
 
 e + e ■- 1.0001 
 
 ,-5 
 
 Starting Vector: x = (0,1) 
 
 Approximate Solution: x = (1.098 x 10" v 9 9.106) 
 
 This example is one indicated by Powell [7] as being difficult, 
 primarily because of being badly scaled. 
 
 Example 8 
 Function: 
 
 f(x) = 
 
 10( 
 
 - x. 
 
 
 ■*2 "I' 
 1 -\ 
 
 Starting Vector: x= (-1.2, 1.0) 
 Solution: x = (l,l) 
 
 This example was used by Powell [7], and appeared earlier in an 
 article by Rosenbrock. [12]. Reference is also made to it by Broyden [ll] , 
 Branin [13], and Fletcher and Powell [8]. 
 
f(x) = 
 
 Ik 
 
 Example 9 
 Function: 
 
 "x 1 (x 1 (5 - x x ) - 2) + x 2 - 13" 
 x 1 (x 1 (l + x 1 ) - lU) + x 2 - 29 
 
 Starting Vector: x = (15, -2) 
 Solution: x = (U,5) 
 
 This problem is given in the Appendix to Powell's article [7] 
 as an example of one which does not converge "by his method. 
 
 Examples 10, 10a, 11, 11a, 12, 12a, 13, 13a 
 
 Function: 
 
 n 
 
 £ (A. . sinx. + B. . cosx.) - E. = 0, i=l,2,. . .n. 
 
 The elements A. . and B. . are calculated by a pseudorandom 
 number generator intended to give uncorrelated numbers from a rectangular 
 
 distribution between -100 and +100. The numbers E. are calculated to 
 
 i 
 
 satisfy the equations for a particular vector x* = (x*, x*, x*,...x*) 
 where each of the components was selected by the pseudorandom number 
 generator from a rectangular distribution between -tt and tt . 
 Starting Vector: x = x* + O.ln, where the components of n = (ru , ru , . . .n ) 
 were selected by the pseudorandom number generator from a rectangular 
 distribution between -tt and tt . 
 
 Solution: x* = ( x* , x*,...x*) as described previously for this functi 
 For Examples 10 and 10a, n=5; for Examples 11 and 11a, n=10; 
 for Examples 12 and 12a, n=20; for Examples 13 and 13a, n=30. The two 
 examples for each value of n have different coefficients obtained by 
 starting the pseudorandom number generator at different points in the 
 
 on, 
 
15 
 
 two cases. Results of the solution of examples of this nature were 
 reported by Powell [7] and by Fletcher and Powell [8]. 
 
 3.2 Discussion of Results 
 
 The results from computer runs of the seventeen examples using 
 the van Me lie method (IN TECH) and the Fletcher-Powell method ( FLPMIN) 
 are summarized in Table I and Table II. Also shown for comparison are 
 the published results from Boggs [H,5] in Table I, and from Powell [7] 
 in Tables I and II for those of the seventeen examples that were used 
 by them. Table I summarizes the results from the examples having just 
 two equations in two unknowns, while Table II includes those examples 
 having from five to thirty equations and unknowns. 
 
 In Table I the convergence criterion for INTECH is of the 
 same form as that used in van Me lie ' s subroutine DIFMF3, and the value 
 is chosen so that the criterion will be more stringent than those 
 (each different) used by Boggs and by Powell in their examples. For 
 FLPMIN the criterion was chosen to be of the same form as that used by 
 Powell since the sum of the squared functions is needed for other 
 purposes in both methods, but the value used was designed to make it 
 about as stringent as the criterion used for INTECH, assuming the error 
 is about equally divided between the two variables. For the examples 
 in Table II subroutines INTECH and FLPMIN were both altered to use 
 exactly the same error criterion as that used by Powell. Such changes 
 appeared to be necessary in order to make valid comparisons for these 
 examples, which involve many more variables than the ones in Table I. 
 
16 
 
 Table I. Comparison of Results for 
 Examples with Two Unknowns 
 
 EXAMPLE 
 NO. 
 
 INTECH 
 
 FLPMIN 
 
 BOGGS 
 
 POWELL 
 
 I 
 
 J 
 
 NFE 
 
 I 
 
 NFE 
 
 I 
 
 NFE 
 
 NFE 
 
 ACC 
 
 1 
 
 16 
 
 1+ 
 
 2k 
 
 35 
 
 105 
 
 
 
 
 
 2 
 
 11+ 
 
 1+ 
 
 22 
 
 2U 
 
 72 
 
 
 
 
 
 3 
 
 21 
 
 U 
 
 29 
 
 27 
 
 81 
 
 17 
 
 36 
 
 
 
 k 
 
 T 
 
 2 
 
 11 
 
 F 
 
 - 
 
 
 
 
 
 5 
 
 16 
 
 1+ 
 
 2U 
 
 26 
 
 78 
 
 18 
 
 36 
 
 
 
 6 
 
 31 
 
 6 
 
 1+3 
 
 F 
 
 • - 
 
 
 
 
 
 7 
 
 ll+6 
 
 26 
 
 198 
 
 F 
 
 - 
 
 
 
 22 3 
 
 10- 10 
 
 8 
 
 12 
 
 3 
 
 18 
 
 66 
 
 198 
 
 
 
 28 
 
 10- 6 
 
 9 
 
 102 
 
 13 
 
 12 8 
 
 30 
 
 90 
 
 
 
 F 
 
 ID" 6 
 
 I - Number of iterations 
 J - Number of Jacobian evaluations 
 NFE - Number of function evaluations or equivalent; a Jacobian or 
 gradient evaluation is counted as n function evaluations 
 n - Number of unknowns 
 
 Convergence Criteria 
 
 INTECH: Z | f X \ < 10 
 i=l 
 
 -6 
 
 i ? -1? 
 
 FLPMIN: Z (f ) < .5 x 10 
 
 i=l 
 
 BOGGS : 
 
 f 1 ! < 10" 5 , i=l,2,...n 
 
 i ? 
 POWELL: Z (f ) < ACC 
 
 i=l 
 
IT 
 
 Table II. Comparison of Results for 
 Examples with More than Two Unknowns 
 
 EXAMPLE 
 
 
 INTECH 
 
 FLPMIN 
 
 POWELL 
 
 NO. 
 
 n 
 
 
 
 
 
 
 
 I 
 
 J 
 
 NFE 
 
 I 
 
 NFE 
 
 NFE 
 
 10 
 
 5 
 
 9 
 
 3 
 
 2k 
 
 33 
 
 198 
 
 11 
 
 10a 
 
 5 
 
 9 
 
 3 
 
 2k 
 
 39 
 
 23U 
 
 12 
 
 11 
 
 10 
 
 13 
 
 3 
 
 1+3 
 
 k6 
 
 506 
 
 19 
 
 11a 
 
 10 
 
 10 
 
 2 
 
 30 
 
 k6 
 
 506 
 
 23 
 
 12 
 
 20 
 
 13 
 
 3 
 
 73 
 
 79 
 
 1659 
 
 36 
 
 12 a 
 
 20 
 
 F 
 
 - 
 
 - 
 
 77 
 
 1617 
 
 36 
 
 13 
 
 30 
 
 F 
 
 - 
 
 - 
 
 107 
 
 3317 
 
 U7 
 
 13a 
 
 30 
 
 19 
 
 3 
 
 109 
 
 113 
 
 3503 
 
 U3 
 
 I - Number of iterations 
 J - Number of Jacobian evaluations 
 NFE - Number of function evaluations or equivalent; a Jacobian or 
 gradient evaluation is counted as n function evaluations 
 n - Number of unknowns 
 
 Convergence Criterion 
 
 For all methods in the table : 
 
 n 
 
 (f 1 ) 2 < 10" 3 
 
 i=l 
 
18 
 
 The subroutine INTECH fails to converge for two of the seventeen 
 cases, Examples 12a and 13 as shown in Table II. It does converge, 
 however, for Examples 12 and 13a, which are of the same nature, involving 
 twenty and thirty variables, respectively, but formed from different 
 sets of random numbers. Apparently the method of determining the values 
 of H and ALPHA to be used is not completely satisfactory in these two 
 cases because the error, after decreasing for awhile, increases markedly 
 for several steps then gradually and continues to cycle in this manner. 
 It was possible to obtain convergence for Example 12a by making a change 
 in the method so that slightly smaller values of H were used. This 
 change increased the number of function evaluations required in other 
 cases, however, and no further attempts were made to optimize the 
 strategy because the primary purpose of the investigation was to 
 evaluate the method as it had been developed up to this time. 
 
 The Fletcher-Powell subroutine, FLPMIN , applied to the 
 function formed by taking the sum of the squares of the original functions, 
 failed to converge to a solution of the original equations for Examples 
 k, 6, and 7- In the case of Example k it converges to a nonzero minimum 
 of the function, which is not a solution of the original equations. In 
 the cases of Examples 6 and 7 it will be necessary for the reader to 
 refer to Fletcher's algorithm [10 ] for the definitions of the names 
 given in quotation marks in the following discussion. In Example 6 the 
 reason for non- convergence was that the variable "gb" as calculated just 
 after the label "search along s:" becomes positive, which means that the 
 matrix "h" is no longer positive definite as it should be if the 
 procedure for updating it were exact [8]. In Example 7 the problem is 
 
19 
 
 caused by continuous looping through the section "beginning with the label 
 "interpolate:", without ever satisfying the criteria for leaving that 
 loop. It is possible that the difficulties encountered in these latter 
 two examples could be circumvented by changes in the subroutine, but no 
 attempt was made to do that in this investigation. 
 
 In the twelve examples in which both INTECH and FLPMIN converge, 
 INTECH requires fewer function evaluations in all except one case. On 
 the average it requires less than one-half of the function evaluations 
 required by FLPMIN in the examples of Table I, which involve only two 
 variables. In the examples of Table II, which involve more than two 
 variables , FLPMIN seems to require on the average about n-times as many 
 function evaluations as INTECH where n is the number of variables in the 
 example. 
 
 Examples 1 and 2 are the same problem with different starting 
 vectors that were expected to lead to different solutions. The same is true 
 of Examples 3 and k. For Examples 1 and 2 INTECH produced the expected solutions 
 in each case while FLPMIN converged to the same solution (the one expected 
 for Example 2) in both cases. For Example 3 it converged to the solution 
 expected for Example h, and for Example U — as shown in Table I — it found 
 a nonzero minimum which is not a solution of the original equations. 
 Actually, it is not surprising that such differences should occur "between 
 the two methods since the function being minimized in FLPMIN is the sum 
 of the squares of the functions used by INTECH, and such a transformation 
 could be expected to produce different regions of convergence to the 
 various solutions. In the practical application of any method to complex 
 nonlinear equations that have not been previously solved, the problem of 
 determining an appropriate starting vector for that particular method 
 will exist. 
 
20 
 
 Results from Example 3 and Example 5 were also given by Boggs [,U,5] 
 Table I indicates that INTECH requires from twenty to thirty percent fewer 
 function evaluations than the Boggs methods. Actually, the table shows 
 the best results obtained by Boggs with each of the two examples out of 
 twelve different implementations of his basic method that were applied to 
 each example. It turns out that his best algorithm for Example 3 did not 
 converge for Example 5 and vice-versa. Considering only those of his 
 algorithms which converged for both of his examples, the best one requires 
 more than twice as many function evaluations as INTECH in each case. 
 
 Referring to the results reported by Powell [T] for Examples 7, 
 8, and 9? it can be seen in Table I that INTECH compares favorably, 
 requiring fewer function evaluations for Example 7 and Example 8, and 
 converging where Powell's method does not for Example 9. Powell's method, 
 however, requires consistently fewer function evaluations than INTECH for 
 Examples 10 through 13a as shown in Table II. Comparisons cannot be 
 completely straightforward in these examples, of course, because the 
 coefficients were generated by pseudorandom number generators and were 
 probably entirely different in the two cases. 
 
 After completing the computer runs on all the examples, it was 
 pointed out to the author that an example from Rosenbrock [12 ] , which 
 was discussed in a recent article by Branin [13], was one for which 
 certain methods did not converge. The reference was found and the example 
 was quickly tried on both INTECH and FLPMIN , requiring twenty-nine function 
 evaluations with the former and 19 8 with the latter. Subsequently, the 
 author realized that this example had previously been done as Example 8, 
 having been found in Powell's article [7] in the first instance. On 
 comparing the number of function evaluations required in the two runs of 
 the same problem, however, it was found that Example 8 only required 
 
21 
 
 eighteen function evaluations with INTECH, although with FLPMIN it 
 required the same number, 198, as on the later computer run. The 
 discrepancy was resolved when it was noted that this example is one in 
 which one of the variables, called x in Example 8, appears only 
 linearly. In doing Example 8 this fact had "been properly noted and the 
 special provision in INTECH for linear variables had been utilized in 
 implementing the example. When the same example was inadvertently 
 used the second time, the linearity of x was overlooked and both 
 
 variables were treated as possibly nonlinear ones in the implementation 
 for INTECH. Thus, the reduction in the number of function evaluations 
 required in van Melle's method when advantage is taken of the special 
 provision for linear variables was documented, although this had not 
 been planned as part of the investigation. 
 
 Example 9, the one for which Powell's method did not converge 
 although success was obtained with both INTECH and FLPMIN, was given as 
 an example in an earlier article by Freudenstein and Roth [l^]. 
 Convergence was obtained by their method, although the total number of 
 function evaluations required was not given. Using a parameter 
 perturbation technique, ten steps were required, each step involving the 
 solution by Newton's method of a nonlinear system having different parameters 
 
22 
 
 k. CONCLUSIONS 
 
 The primary conclusion is that the integration technique 
 developed by van Melle for the solution of nonlinear equations has 
 considerable merit. As implemented here in the subroutine INTECH , it 
 was tried on a number of different examples. Most of these had been 
 identified in previous publications as examples which were known to 
 be troublesome to solve, at least by some methods. Every example that 
 was tried on INTECH is reported in Section 3. Under these conditions, 
 the fact that convergence was obtained in fifteen of the seventeen 
 cases tried, and that the number of function evaluations was less than 
 for the other methods compared in about half the cases , gives support 
 for this conclusion. 
 
 Another conclusion, based on thinking about the situation 
 rather than on experimental data, is that it is probably feasible to 
 further improve the implementation of this method. The procedure used 
 to establish the step size H and the parameter ALPHA for the next 
 iteration was determined empirically; and since it works well in many 
 cases , it should be retained in basically the same form unless 
 extensive additional testing proves some different approach to be better. 
 It is suggested that by studying in detail all examples for which the 
 method is found not to converge, or to converge slowly, it may be 
 possible to supply additional logic to recognize in the early stages 
 similar cases and take more intelligent action for each of such 
 specific situations studied. Examples in which these particularly 
 troublesome situations did not develop would be handled by the 
 present procedure. 
 
23 
 
 The empirical procedure for deciding when to re-evaluate the 
 Jacobian and its inverse seems to he reasonably effective in reducing 
 the number of function evaluations in INTECH , as indicated particularly 
 by comparing with the results from FLPMIN in Table II, which includes 
 all the examples involving five or more variables. However, it is 
 possible that matrix updating procedures such as those described by 
 Broyden [15] or the one used by Powell [6,7] would prove to be even 
 more economical of computer time, especially for systems with many 
 unknowns. This possibility is supported by the fact that in the 
 examples of Table II, INTECH consistently requires more function 
 evaluations than Powell's method, and the number of these used for 
 the Jacobian was greater than the number used for iterations. It is, 
 therefore, concluded that the use of matrix updating procedures in 
 combination with the basic approach of van Melle may be well 
 worth investigating. 
 
2k 
 
 LIST OF REFERENCES 
 
 [l] van Me lie, Bill, "AN OPERATING MANUAL FOR DIFMF3: A Program to 
 Solve Initial Conditions in Simultaneous Differential 
 Equations ," Department of Computer Science File UIUCDCS-F-72-870 , 
 University of Illinois at Urb ana-Champaign , April 1972. 
 
 [2] Gear, C. W. , "Simultaneous Numerical Solution of Differentials- 
 Algebraic Equations," IEEE Transactions on Circuit Theory , 18 
 #1, 1971, PP. 89-95. 
 
 [3] Gear, C. W. , NUMERICAL INITIAL VALUE PROBLEMS IN ORDINARY 
 
 DIFFERENTIAL EQUATIONS, Prentice-Hall, Inc.: New Jersey, 
 1971, PP. 223-227. 
 
 [h] Boggs , P. T. , "The Solution of Nonlinear Systems of Equations by 
 A-Stable Integration Techniques," SIAM J. Numer. Anal. , _8, 
 #U, 1971, PP. 767-785. 
 
 [5] Boggs, P. T. , "The Solution of Nonlinear Operator Equations by 
 A-Stable Integration Techniques," Doctoral Thesis, Cornell 
 University, Ithaca, New York, 1970. 
 
 [6] Powell, M. J. D. , "A Hybrid Method for Nonlinear Equations," in 
 NUMERICAL METHODS FOR NONLINEAR ALGEBRAIC EQUATIONS, ed. 
 P. Rabinowitz, Gordon and Breach Science Publishers: 
 New York, 1970, pp. 87-ll 1 4. 
 
 [7] Powell, M. J. D. , "A FORTRAN Subroutine for Solving Systems of 
 Nonlinear Algebraic Equations," in NUMERICAL METHODS FOR 
 NONLINEAR ALGEBRAIC EQUATIONS, ed. P. Rabinowitz, Gordon 
 and Breach Science Publishers: New York, 1970, pp. 11 5-1^-9 • 
 
 [8] Fletcher, R. and Powell, M. J. D. , "A Rapidly Convergent Descent Method 
 for Minimization," Comp. J., 6, 196 3, pp. 163-I68. 
 
 [9] Wells, M. , "Algorithm 251 Function Minimization [EU]," Comm. ACM , 
 8, 1965, pp. 169-170. 
 
 [10] Fletcher, R. , "Certification of Algorithm 2 51 [EU] Function 
 Minimization," Comm. ACM , 9, 1966 , pp. 686-687. 
 
 [ll] Broyden, C. G. , "A New Method of Solving Nonlinear Simultaneous 
 Equations," Comp. J. , 12, 1969, pp. 9^-99- 
 
 [12] Rosenbrock, H. H. , "An Automatic Method for Finding the Greatest 
 
 or Least Value of a Function," Comp . J . , 3, I960 , pp. 175-1 8U. 
 
 [13] Branin, F. H. , Jr., "Widely Convergent Method for Finding Multiple 
 Solutions of Simultaneous Nonlinear Equations," IBM J. Res . 
 Develop . , 16 , #5, September 1972, pp. 50U-522. 
 
25 
 
 [ik] Freudenstein, F. and Roth, B. , "Numerical Solutions of Systems 
 of Nonlinear Equations," J. ACM , 10 , 1963, pp. 550-556. 
 
 [15] Broyden, C. G. , "Recent Developments in Solving Nonlinear 
 
 Algebraic Systems," in NUMERICAL METHODS FOR NONLINEAR 
 ALGEBRAIC EQUATIONS, ed. P. Rabinowitz, Gordon and Breach 
 Science Publishers: New York, 1970, pp. 61-73 . 
 
26 
 
 APPENDIX 
 
27 
 
 ) , YLSV<2) ,OEL 
 
 M * EXAM°LE g 
 
 PFAL*S Y(2),nY(2),Fl( 2),YL(2) ,SAVE(2, 
 PEAL*3 YD<2) ,PWQRK(2) 
 
 DIMENSION PW (?,?) 
 N C X = 8 
 
 Y<1 ) =-1.2 DO 
 YL( I )=1.00 
 M=2 
 NL = 1 
 
 DEL=1.0-6 
 NIY = N-NL 
 NSFMD=110 
 WRITE (6,1 J MEX 
 1 FOPMAT (i-PX&vp|_c Nn.«,i c ///) 
 
 CALL INJTFCH{Y,YL,nY,PW,DPL,Pl,YO, S A VE , YL S V , P WOR K , % 
 
 NtNYfNL tNSEND) 
 END 
 SURROt) T I NF HI F«=UN ( D Y, Y , YL , N, NY t NL ) 
 PEAL*3 DY(M), Y(NY) ,YL(NL) 
 
 DY(1 )=10.D0*{ YL(1 )-Y<l )**2) 
 DY( 2)=1.D0-Y( 1) 
 Q PTLPN 
 END 
 SUBROUTINE MAT SCT{ PW,Y, YL t N,NY,NL) 
 REAI *=> Y(MY) , YL ( NL ) 
 DI MENS ION PW(N,N) 
 
 PW( 1 ,1 ) =-20.D0*Y(l ) 
 pwi ! ,2) = !0.no 
 
 PW( 2, 1 )=-l .00 
 PW( 2,2) =0.00 
 
 R E T UP N 
 
 END 
 
28 
 
 SUBPDU T INF INTECHi YtYUDY, PWtDEL t p l tYD,SAVE, YLSV% 
 , PWO°K,N»NY, NL.NSEND) 
 
 REAL*B DY(N),F1(N),Y(NY),YL(NL),YD(NY),SAVE(2,NY),YLSV(NL)?; 
 , PrfORM N) ,0MAX1,DMIN1 , DEL, ALPHA, H f R,DET,SStS ,HOLD ,DD , D , HH 
 
 DIMENSION PW(N,N) 
 
 11 ************************************************************ 
 •« PARAMETFR LIS T 
 
 " Y V=rTRR OF M1NLINEAP VARIABLES (INPUT). 
 
 » YL VECTOR OF LINEAR VARIABLES (INPIJT). 
 
 ii ny VFC T 0» O c FUNCTION VALUES, CALCULATED 
 
 " RY TH^ SUBROUTINE DIF^UN. 
 
 " PW JACOB I AN MATRIX CALCULATED BY SUBROUTINE MATSET 
 
 " ns th c INVERSE JACORIAN AFTER CALLING M[NV. 
 
 » DFL FRROP C.p ITE" IOM.SUOPLIEO BY CALLING PROGRAM. 
 
 •» F! VECTOR EDUAL TO PRODUCT OF PW AND DY . 
 
 ii yd TH C VECT-^P D^ CHANGES TQ BE APPLIED TO VECTOR Y. 
 
 " SAVF TWD-POW VECTDD USED TO ST^PF VALUES OF Y AND YD 
 m rc^ORE niexT IT c 5^j^ f j, i nj C4s c RF-S T ART N r = DED. 
 
 " YLSV V = CTnc? USFD TO STOPE VALUES O c YL BE c OPP 
 " NEX T ITERATION, IN CASE PE-STAPT NEEDED. 
 
 •• PwnPk USED "INLY c OR STORAGE S p ACE WITHIN SUBROUTINE minv. 
 
 * N TUP TDTAL NUMBER OF VARIABLES (INPUT). 
 
 " NY T" c NUMBER OF N1NLIN P AP VA°IABLES (INPUT). 
 
 " NL TUP NUMB C R O c LINEAR VARIABLES (INPUT). 
 
 '• NS p ND THE maximum NUMBER DF ITERATIONS DESIRED (INPUT). 
 
 !•*****************************#-****************************** 
 l» ********** 
 
 " INITIALISE VARIABLES 
 11 
 
 \JD C AT! =1 
 JACOR=0 
 
 p =1 
 
 N 7 - M * N 
 NJ=NY*2 
 
29 
 
 I F(MJ.L T .!0) NJ = 10 
 
 NJJ=NJ*2/3 
 
 ALPHA=0 
 
 H=l 
 
 MS=0 
 
 NW=0 
 i» ******* *** 
 
 '» PQINT HEADINGS 
 
 M 
 
 MRITE(6,1) N,NL,OEL 
 
 1 FORMAT(»iN =»,I4,« NL =»,I^,« DEL ='»D10.2) 
 WRI TE(6,2) 
 
 2 FOR^ATCO MS NW ALPHA H ', 8X, • ERROR ' t 6X, % 
 
 » Y( J) AND YL(* ) •// ) 
 
 II ** ** ****** 
 
 " INITIAL EVALUATION n F FUNCTION (DY). 
 if 
 
 NS = NS+! 
 
 CALL DIFRJN( DY,Y,YL,N, NY, ML ) 
 ii ********** 
 
 " INITIAL EVALUATION qf JAC06IAN (PW). 
 ii 
 
 CALL MATSFT(PW,Y t YL, N f NY, NL I 
 
 ************ 
 
 »• INITIAL c VAH)ATinN| If INV^S^ JACO^IAN ( PW ) . 
 ii 
 
 NW=NW+l 
 
 CALL MTNV( {J W» ! DET,N,NZ,Fl,PWORK) 
 
 « * * ******** 
 
 " PRODUCT °W*ny IN V^CTPR Fl. 
 •• 
 
 C A L 1 M\TMUL(PW,DY»F1, N J 
 •i ******* * ** 
 
 " NORM STORED IN SS FOP STATING VALUES. 
 
 03 ? J = ! t N 
 3 SS=SS+DARS ( ^Y ( J > ) 
 ii********** 
 
 •• pcinjt RFSULTS OBTAINED 4ITH STARTING VECTDR. 
 ii 
 
 WP IT C ( s,i ? ) NS,NW, ALPHA,H,SS,(Y( J) , J=1,NY) 
 I c (Nt. .GT.fl) WRITE! 6 t 1 -^) ( YL( J) , J = 1,ML) 
 
 II********** 
 
 ii PREPARE FOR THF INI T IAL PREDICTOR STEP FOR Y; DO THE 
 
 " INI T IAL CORRECTOR STFP fhr YL . 
 
 •i 
 
 DO 4 J=! , NY 
 
 U Y^( J )=-M*ci ( j ) 
 
 T c < N L . L r . ) GO TO 6 
 
30 
 
 oo 5 J=1,NL 
 
 5 YL( J)=YL( J)-F1 ( J+NY) 
 
 6 CONTINUE 
 
 ********** 
 
 SAVE PRESENT VECTORS Y AND YL. 
 
 DP 7 J=1,NY 
 SAVEd , J» = Y( J ) 
 
 7 SAV C (2, J)=YD( J) 
 I C (NL .L r .O) G^ TO 9 
 o fl J = ! f N L 
 
 R YLSVf J) =YL( J) 
 
 9 CONTINUE 
 
 ********** 
 
 INOIC4 T F NEWTON'S METHOD IS IN USE BY NEWTON = 1. 
 
 NEWTQN=! 
 ********** 
 
 BEGIN MAIN ITERATION LOOP BY COMPLETING THE PREDICTION OF Y. 
 
 10 DO I 1 J = l ,NY 
 
 11 Y( J)=Y( J)+YD< J) 
 
 ** ******** 
 
 FVALUATE T HF FUNCTION (DY). 
 
 N S = N S + 1 
 
 CALL DIP-UN(DY,Y f YL,N T NY ,NL ) 
 ********** 
 
 CALCULATE THE VALUE OF T HP N n <5M (S). 
 
 S=0 
 DO 1 -> 
 
 uu J = 1 ,N 
 
 12 S=S*0ABS<DY( JH 
 
 ********** 
 
 ********* 
 
 PRINT RESULTS FOLLOWING PREDICTED STEP, THEN RETURN IF 
 NORM IS L C SS THAN QR c OUAL T D C L. 
 
 WRI T >-(^, !3 ) NS,NW, ALPHA, H, S,(Y(Jl,J=l,NY) 
 13 FORMAT (14, I3 v F6«2«2fUl«2t2Dl5.5/(5X v 4D15«5)) 
 
 IF(NL.ST.O) WRITE! 6, 14) ( YL( J) , J--!,NL) 
 ] U F 00 MAT (^X,^D!5.5) 
 
 IF(S.LF.DFL) RETURN 
 
 ** ******** 
 
 CALCULATE PREP (PATIO OF NO p M/ °R c V I OUS NORM). 
 
 PR ED=S/SS 
 ********** 
 
 BEGIN PROCESS OF DETERMINING VALUES op ALPHA AND P 
 BY WHICH H WILL BE M UL.TIPLIEO) e ?R N-XT I T c RATION. 
 
 ( FACTOR 
 
31 
 
 GO 
 
 TO 
 
 35 
 
 on 
 
 T l 
 
 IS 
 
 GO 
 
 TO 
 
 30 
 
 IF ( NEWTON. EO. 1) 
 IF< PREO.LT.100) 
 IF(NOFAIL.Eg.l ) 
 
 !5 CONTINUE 
 
 IF (PPE0.LT..9R) GO TO 16 
 
 ALPHA=l 
 
 PaDMiMl ( 1. 300, 0. 600/H ) 
 
 GO to 17 
 
 16 D =1.7DT-.e5D0*H+,! 5D0/H 
 ALPHA=ALPHA*0.8 
 
 " SAVF PRESENT VECTORS AND VALUE 
 it 
 
 17 DO 19 J=1,NY 
 <^AVE( 1, J)=Y( J) 
 
 IB SAVE(?, J)=YD( J) 
 
 IF(NL.LF.O) GO T 2 
 
 0^ 19 J=1,NL 
 19 YLSV( J>=YL( J) 
 2 CONTINUE 
 
 HOLD=H 
 
 •« ********** 
 
 ri D c T cpMlNr W HETHEP J AC OR IAN SHOULD 
 IF( SS.GF.l.OOJ GO TQ 21 
 !»=( JACTR.LT.NJJ ) 00 T3 ?-> 
 21 CONTINUE 
 
 JAC 3B = JAC0R-1 
 
 IF( JAC">B.G T .0) GO TH ?3 
 
 ************ 
 
 " EVALUATE TH C JACOBIAN (OW), 
 
 OF H BEFORE NEXT ITERATION. 
 
 BE RE-EVALUATEO. 
 
 22 CALL MA T S- T (PW,Y,YL,N,NY,NL) 
 
 M w=N W+l 
 
 »♦ * * * * * * * * A * 
 
 " EVALUATE the INVERSE QF JACORIAN 
 
 (PW) . 
 
 CALL MINV(PWfOET,N t N/,Fl,PwnPK) 
 J A C B = N J 
 
 f D m v e c t o p 
 
 THE PonrjuCT PW*OYS THEN DO THE 
 
 CORRECT^ S Tr P ^OP Y. 
 
 23 CALL M* T M|jL( ^,nY, c l, N) 
 IF( ALPHA. LT. 0.01 ) GO Tn 
 DD=1 .D0/( ! *H*ALPHA > 
 qo -»4 J = l , NY 
 n=( c l ( I)*H*YD(J) )*D0 
 Y( I )=Y( J}- ALPHA*D 
 
 ?4 YD{ J)=P*(YO(J )-Q) 
 
 GO t n -> 7 
 
32 
 
 25 HH=-H*R 
 
 nn 26 J=l ,NY 
 
 26 Y0< J)=HH*FK J) 
 
 27 CONTINU c 
 IP(NL.LE.O) GO 
 
 " DO T HE CORRECTOR 
 
 TO 2o 
 
 STEP FDR THF LINEAR VARIABLES IN YL. 
 
 on 2 8 J = 1 ,NL 
 2B YL( J) = YL( J )-Fl (J+NY) 
 
 29 SS=S 
 
 NOFAIL-1 
 
 I c ( NS.GE.NSENO) 
 
 " END D c mmm 
 
 GO TO 3 7 
 nation LOOP. 
 
 GO TO 10 
 30 CONTINUE 
 
 *• CHANGE AL p HA f R, 
 
 ALPHA=! .0 
 IP(JACOB.L T .NJJ 
 Q = r>M AVI ( 0.500, 
 
 H=HOLD*P 
 
 " PREVIOUS VFfTPRS 
 " RE-STARTING LAST 
 
 OP 3 2 J=l ,NY 
 Y< J)=SAVE(1 , J) 
 Y0( J )=9AVF( 2, J) 
 IF(NL-L c .O) GO 
 on ^3 |=t,NL 
 YL( J)=YLSV< J ) 
 FONT I NUE 
 MOFAIL=0 
 GO TO ! 
 CONTINUE 
 
 If********** 
 
 " IF ( N F W TH N s ] A N n 
 
 " C0PR C CTI0N PROCES 
 
 IF(P»E0,LT.0.95 
 
 ii *# *^ * ** a* * 
 
 " BUT IP PRED>=.9 5 
 " NEWTON, J AC 01; TH 
 '• ALTFP VALUE n F NO 
 
 *1 
 
 3-> 
 
 ^3 
 3^ 
 
 •*«; 
 
 ANO H C 0R THE CASE (PRFDMOO AND NOFAIL=l) 
 
 ) JACOB=0 
 
 . 7 D0/HOL n ) 
 
 RESTORED AND NOFAIL S r T TO ZERO BPF3RF 
 STcp, BECAUSE n F LARG- INCREASE IN ERROR 
 
 T 3^ 
 
 PRFD<. : >E) GO BACK TO 00 THF NORMAL 
 S, REMAINING IN THE NEWTON METHOD. 
 ) GO TO 17 
 
 PRIMT A MFSSAG-, CHANGE ALPHA, H, R, 
 EN GO TO tFSTORE PREVIOUS VECTORS AND 
 FAIL BEFORE RE-STARTjsjG LAST ITERATION, 
 
 WRI T=(6, ->6)S, ( Y< Jl , J = l ,NY) 
 
33 
 
 3 6 poo 
 
 IH 
 
 H=. 
 
 ALP 
 R=. 
 JAC 
 MP W 
 GO 
 37 COM 
 n ****** 
 
 " IF TH 
 
 " PRTMT 
 
 W&I 
 
 3 8 
 39 
 40 
 
 FOR 
 WRI 
 POR 
 WRI 
 FOR 
 RET 
 EMO 
 
 MAT ( 
 
 NL.OT 
 
 01 
 
 HA = 1 
 01 
 0R = 
 
 TOM = 
 
 m 3i 
 
 TIMUF 
 **** 
 
 F MUM 
 
 THP 
 TE(6, 
 MAT( • 
 TE<5, 
 
 M \ T ( • 
 
 TE(6, 
 MAT ( 
 
 UPM 
 
 • NEWTON FAILED: • ,9X,011.2, 201 5 . 5 / < 5X ,4015 . 5 ) ) 
 .0) WRITE(6,14) (YL(J) ,J=1 ,NL) 
 
 BE c? nc FUMC T IOM EVALUATIONS HAS 
 
 VEC T ORS DY, Y f YL , !\N0 F!' t THEN 
 
 3«) (0Y{ J) ,J=1 ,N) 
 
 -DY: •/( •**• ,4014.6) ) 
 
 ^ Q ) I Y( J), J=T ,NY) 
 
 -Y: •/(•** ' ,4016. 9. ) ) 
 
 40) ( J,F1( J),J=1,N) 
 
 »-F! :«/31 I5,«**« ,014.6) ) 
 
 EXCEEDEO NSFND, 
 3FTIJ9N. 
 
 M ***** * 
 H ** **** 
 •f 
 
 SUBP0UT 
 
 REAL+3 
 01 MEMS 
 
 no 
 
 SljM 
 
 on 
 
 sum 
 
 FT ( 
 RF T 
 
 en 
 
 40 
 50 
 
 ************************* ********* ********************* 
 * * ***************************************************** 
 
 rfsJE MATMUL(PW,DY , c l ,N) 
 
 0Y( M) ,F1( M) , SUM 
 ION! PW(N,N) 
 *0 J-l ,M 
 = 
 
 40 K=l , M 
 
 = SIIM+PW( J ,K) *0Y( K) 
 J ) = S!JM 
 
 URN 
 
FormAEC-427 U.S. ATOMIC ENERGY COMMISSION 
 
 aecm3201 UNIVERSITY-TYPE CONTRACTOR'S RECOMMENDATION FOR 
 
 DISPOSITION OF SCIENTIFIC AND TECHNICAL DOCUMENT 
 
 ( See Instructions on Reverse Side ) 
 
 1. AEC REPORT NO. 
 
 C00-li+69-0226 
 
 2. TITLE 
 
 EVALUATION OF AW INTEGRATION TECHNIQUE FOR THE 
 SOLUTION OF NONLINEAR EQUATIONS 
 
 3. TYPE OF DOCUMENT (Check one): 
 
 QJ a- Scientific and technical report 
 
 I I b. Conference paper not to be published in a journal: 
 
 Title of conference 
 
 Date of conference 
 
 Exact location of conference. 
 
 Sponsoring organization 
 
 □ c. Other (Specify) 
 
 4. RECOMMENDED ANNOUNCEMENT AND DISTRIBUTION (Check one): 
 
 [3g a. AEC's normal announcement and distribution procedures may be followed. 
 
 | | b. Make available only within AEC and to AEC contractors and other U.S. Government agencies and their contractors. 
 
 i J c. Make no announcement or distrubution. 
 
 5. REASON FOR RECOMMENDED RESTRICTIONS: 
 
 6. SUBMITTED BY: NAME AND POSITION (Please print or type) 
 
 C. W. Gear 
 
 Professor and Principal Investigator 
 
 Organization 
 
 Department of Computer Science 
 University of Illinois 
 Urbana, Illinois 6l801 
 
 Signature 
 
 ^vfil^^-. 
 
 Date 
 
 August 19 73 
 
 FOR AEC USE ONLY 
 
 7. AEC CONTRACT ADMINISTRATOR'S COMMENTS, IF ANY, ON ABOVE ANNOUNCEMENT AND DISTRIBUTION 
 RECOMMENDATION: 
 
 8. PATENT CLEARANCE: 
 
 I I a. AEC patent clearance has been granted by responsible AEC patent group. 
 [~~l b. Report has been sent to responsible AEC patent group for clearance. 
 I I c. Patent clearance not required. 
 
JIBLIOGRAPHIC DATA 
 HEET 
 
 1. Report No. 
 
 UIUCDCS-R-73-576 
 
 , Title and Subtitle 
 
 EVALUATION OF AN INTEGRATION TECHNIQUE FOR THE 
 SOLUTION OF NONLINEAR EQUATIONS 
 
 3. Recipient's Accession No. 
 
 5- Report Date 
 
 August 1313 
 
 Author(s) 
 
 James Harold Robb 
 
 8- Performing Organization Rept. 
 
 No. 
 
 Performing Organization Name and Address 
 
 I Department of Computer Science 
 
 University of Illinois 
 Urbana, Illinois 618OI 
 
 10, Project/Task/Work Unit No. 
 
 11. Contract/Grant No. 
 
 US AEC AT(ll-l)lU69 
 
 2. Sponsoring Organization Name and Address 
 
 US AEC Chicago Operations Office 
 9800 South Cass Avenue 
 Argonne, Illinois 
 
 13. Type of Report & Period 
 Covered 
 
 Thesis Research 
 
 14. 
 
 5. Supplementary Notes 
 
 6. Abstracts 
 
 The method for solving systems of nonlinear algebraic equations 
 that is the subject of this investigation was developed by van Me lie [l]. 
 His subroutine DIFMF3 was designed to solve the steady-state problem of 
 ordinary differential equations, and for this purpose each component of 
 the vector to be solved for could be either a derivative or an ordinary 
 variable. The development of this method grew out of ideas presented by 
 Gear [2,3] concerning the solution of systems of mixed algebraic and 
 differential equations. 
 
 7. Key Words and Document Analysis. 17a. Descriptors 
 
 nonlinear 
 algebraic equations 
 
 7b. Identif iers/Opcn-Ended Terms 
 
 |7e. COSAT1 Field/Group 
 
 B. Availability Statement 
 
 unlimited 
 
 19. Security Class (This 
 Report) 
 
 UNCLASSIFIED 
 
 20. Security Class (This 
 
 Page 
 UNCLASSIFIED 
 
 21. No. of Pages 
 
 39 
 
 22. Price 
 
 3RM N TIS-'JS (10-70) 
 
 USCOMM-DC 40329-P7 1