Computer realization of the method of Newton diagrams

6 REFERENCES 1. KANTOROVICH L.V., Functional analysis and applied mathematics, Usp. Mat. Nauk, 3, No.6, 89-185, 1948. 2. LUCHKA A.YD., Acceleration of the convergence of gradient methods, in: Mathematical software and efficient organization of the computational process, No.2, In-t Kibernetiki Akad. Nauk UkSSR, Kiev, 58-65, 1969. (Priblizhennoe 3. KRASNOSEL'SKII M.A., et al., Approximate solution of operator equations reshenie operatornykh uravnenii), Nauka, Moscow, 1969. 4. LUCHKA A.YU., Projection-iterative methods of solving differential and integral equations (Proektsionno-iterativnye metody resheniya differentsial'nykh i integral'nykh uravnenii), Naukova Dumka, Kiev, 1980. 5. KANTOROVICH L.V. and AKILOV G.P., Functional analysis is normed spaces (Funktsional'nyi analiz v normirovannykh prostranstvakh), Fizmatgiz, Moscow, 1959. 6. EPMAKOV V.V., and KALITKIN N.N., Two-stage gradient descent, Zh. vych. Mat. i mat. Fir., 20, No.4, 1040-1045, 1980. 7. SAMAPSKII A.A. and NIKOLAEV E.S., Methods of solving mesh equations (Metody resheniya setochnykh uravnenii), Nauka, Moscow, 1978.

Translated

U.S.S.R. Comput.Maths.Math.Phys Printed in Great Britain

by D.E.B.

0041-5553/84 $lO.OO+O.oO 01985 Pergamon Press Ltd.

.,Vol.24,No.4,pp.6-14,1984

COMPUTER REALIZATION OF THE METHOD OF NEWTON DIAGRAMS* A.M.

BALK

A FORTRAN program is given for finding implicit functions as generalized power series. A generalization of the method of Newton diagrams is described for solving differential equations.

1. Let The equation

f(z,w)

be an analytic

function

of two

(complex or real) variables,

f(0,O)=O. (1)

l(z,w)=O

in general defines in the neighbourhood of the point z=O several continuous implicit functions w(z) such that u)(O)=O. In practice we are often faced with the need to find these functions, e.g., in connection with the bifurcations of stationary solutions of differential equations, see /l/. Newton gave a method whereby all the solutions w(z) can be found as series in fractional powers of z: lu(z)=Rz'+R'z"+R"z""+ .... (2) where e, E',~~,...,is an increasing sequence of positive rational numbers. Note that the more general case, when w is an n-dimensional vector, and f is an n-dimensional vector function, can often be reduced to solving the scalar equation (l), see /2, 3/. Below we given a FORTRAN program realizing the method of Newton diagrams. 2.

Let us first recall the essentials

of Newton's

method.

where all a,(z),n=O, 1,.... are analytic functions. We write identically, as a,(z)=a,~'~+a.,z'"'+. . . . where all coefficients Q.~ are non-zero, T..r
while r*e,.. .,

the powers

We write

f(z,~) as

....

f(z,~~)=s~(z)+a,(z)w+ . ..+a.(z)w”+

the

(3) a.(Z) which do not vanish

r-k are non-negative

+f”D... .

there is necessarily one or more (or an infinity of) least numbers; In other words, we choose m such that least, call it m.

integers,

and (4)

from these we choose the

(5) rro>ymo VnO, RM, W,==w,(z)+O, z+O. and R from the condition that, on substituting the solution w(z) given by (6) in the function After this substitution, f(z,u)) takes instead ofw, we must obtain identical zero. f(z,w), the form

lZh.vychisl.Mat.mat.Fiz.,24,7,972-985,1984

7

f (z, w) = (a,&“fa*,z~“+

. . .) +(a,&‘~*fa,,z’~+

+ . . . i(Q,,Z”*i-amrz’“~+.

. . .) [z’(Rfw,)]

1n particular, the sum of "1owest"termsmustvanishin lowest order of smallness as z-&). For any fixed e>O certainly be found among the following:

. .) [z’(R+luA I”+ . . . .

thisexPression (i.e., thetermshaving all the lowest terms of (7) will

(8)

C&z", a,&*(Rz'),....amoz'~(Rz')"'. The powers

(71

to which z is found in these terms are yao,rlo+E,...,

(9)

rm,+ne.

Since the sum of all the lowest terms in (7) must vanish identically, at least two of powers (9) must be equal, while the rest must be not less than these. In other words, e has to be found from the conditions: there exist two numbers i, j.i#j, OCi,jO. Newton gave a geometric method of finding E from this non-singular system of equations and inequalities. For the numbers n, n
(11)

h rn.).

The set of all these Points is called the Newton diagram. With every point of (11) we .asso-ciateanumber omr (see Fig.1). The absence on the axis of ordinates of a point of this diagram obviously implies that c(z)-0 , and hence the function (Leo is a solution of (l), while the multiplicity of this solution is equal to the least n for which a point (11) is defined. To find the non-zero solutions, we can use the following arguments. We fix a number e>O. Through all the points of the diagram we draw straight lines with slope e (we understand this as follows: the tangent of the angle between any such straight line and the negative direction of the abscissa axis is e, see /4/ and Fig-l). The line, passing through the point Hence the (n, Y=.). cuts the ordinate axis at a height h.- r.,+ne (above the abscissa axis). problem of solving system (10) can be restated as: to find an e>O such that, fortwodistinct numbers i, j, OCf, j
n

L?

Ill

Fig.1

Fig.2

We construct the convex hull of all points of the diagram (Fig.2). Let L be its "southwest boundary" (i.e., the set of segments of the boundary with the following typical property: a straight line containing the segment cuts both sides of the first coordinate angle and separates the diagram fran the origin). Then every link L,of step-line L defines the required We call L power s: it is equal to the slope of link L,. Obviously, e is a rational number. the Newton step-line. We now proceed to find the coefficient R in expansion (6). Let link L, contain r points (r*2) :

(h, -fd, . . . , (i., rd. This means that the powers are equal to one another, while the lowest terms in (8) are

y&ire,..., r,,+i,e all the remaining powers

in (9) are greater

than them.

Thus

cr(,P"(R~~",...,~z'h((~*~t Their sum must vanish; since braic equation in R,%%:

z appears in these terms

to the same power, we obtain

(12)

aJz~+...+a&-o. On solving this equation,

we find all the possible

the alge-

values

of R corresponding

to the link

L,. Hence, on inspecting all the links of step-line L, we find all the possible pairs (e,R) Every such pair defines the first (to certain c there may correspond several distinct R). term in the expansion (2) of a solution of (1). To find the next term in (2), it is best to use the change of variables z-z,Q, W-Z$P(R+W,),

(13)

8 where P/Q is the irreducible fraction equal to a; for variables z, and&we obtain of the same type as (1). tion f,(Z,, Wi)-0 The following propositions, proved in the past by various authors, see /2-4/, tial for justifying and applying Newton's method. Proposition

1.

Proposition Newton's method,

Equation

(1) has only a finite number of continuous

the equaare essen-

solutions.

2. There is an integer Q0 such that, in all expansions (2) obtained by the powers E, E',... can be written as fractions with denominator Qo.

Proposition 3. The sum of the degrees of the algebraic equations, from which all poSsible values of R' can be found (see (2)), does not exceed the multiplicity of the root R. Hence, if R is a simple root, the pair (e',R') is uniquely defined,fl being given by a linear equation. Proposition 4. Newton's method gives all the solutions of (11 in a neighbourhood of the point z=O. Proposition 5. All the series obtained by Newton's method are convergent in a neighbourhood of the point z=O and solutions of (11.

3. series,

A computer realization i.e., polynomials.

12345

of Mewton's

method works only with finite segments

of a power

It is important to choose a reasonably good representation of the polynomial in the computer memory. It is unsatisfactory to use the natural matrix form in which the term z'w' of the polynomial appears at the (i,i)-th position, since, when finding the next term of expansion (2), we have to make substitution (13). If the polynomial f(z,w) contains the monomial a+?'w', then the polynomial f,(z,,mI) obtained from f(z,w) after substitution (13), will in

2 particular contain the term CZ,,Z,*+'~UJ?, Hence it can be seen that the dimensionality of the matrix of coefficients a# a,(z)au(z) will increase rapidly (when finding successively the terms q@,@(z)-0 of the expansion of the function u)(z)), while the matrix itself will contain many zeros. We thus preferred a Fig.3 different form of the polynomial for the computer realization of Newton's method. Two large one-dimensional blocks of the same size are introduced: the integral block 2 and the real (or complex) block A. Also, an integral one-dimensional block of powers E'L (from the French word fleche = arrow) was introduced, the length of which is 1 greater than the degree of the polynomial f(z,m) in zu (Fiq.3). Notice that substitution (13) does not increase the degree of the polynomial in u) (to be more precise, the degree of the polynomial f(z,W) in w). We start by registering succesfrh, l4) in w,is not higher than the degree of sively in blocks Z and A the information about each term of the polynomial which contains W to the zeroth degree (the order of registering these terms is unimportant): Z(l) (i.e., the content of cell J of block Z) is the degree of this term in z, and A(I) is the coefficient of The number of the cell in which the last term with zero degree in UJ is written this term. is registered in the first cell of block FL; if the polynomial has no terms containing W to (in other words, in the first cell of block FL is indicated the zero degree, then FL(l):=0 Thus, in number of all non-zero terms of the polynomial f(z,W), not containing mat all). the cells of blocks Z and A with numbers from 1 to FL(l) is written the polynomial G(Z), see (3). Then, we register in blocks Z and A the information about the non-zero terms, containing (D to the first degree, i.e., the information about the polynomial a,(z) : the number of the cell containing the last term of this polynomial is written in the second cell of block If, for some k, there are not terms at all in the polynomial f(z,w), that contain FL, etc. (k+l)-th) cell of block FL we register the content of the prethen in the next (i.e., Wk, vious cell of this block (if there are no terms in the polynomial whose degree in wis less than some number s, then FL(s):=...=FL(i):=O). Thus, the complete information about polynomial f(z,w) is written in the three blocks Z, A, FL. Let us briefly describe the program realizing Newton's method. A

Algorithm

1.

At the program input the three blocks DA, stage 1: input and transformation of data. DZ, DW (where D denotes "data") are transmitted, containing information about each monomial of the initial polynomial, namely: in the cell of these blocks having the same number are located respectively the coefficient of the monomial, its degree in z, and its degree in W (we assume that like terms are collected in the polynomial and that terms with zero coeffiOperation of the program cients are not written). The order of the monomials is arbitrary. starts with finding the degree of the polynomial in LU, while the polynomial itself is located in blocks 2, A, FL. Running over the elements of block 2 in stage 2: construction of Newton's diagram. turn, we find the order-wise first least number in it; this is the number which was called in sect.2 (knowing the number of the cell in which it is written, we can use block FL rm The cell number of block 2, to find the number m itself; but this will not be required). in which yno is found, is written in a special cell, call it HL. We reduce the content of each cell of block Z by yrnp: this implies division of /(~.a) by zln*.

9

In Newton's

method

a special

role is played by the monomials,

see Sect.2, a.~zTrw", ?&in. n,(z)+O,) On sorting successively through the elements of block 2 with numbers not exceeding HL, we find the triples n
(We recall that the number n only runs over the values for which

solution Stage 3: study of Newton's diagram. Let us see if (1) has the zero w(z)-0. If X(1)-0, then Newton's diagram includes the point (O,yoO), lying on the ordinate axis, and hence there are no zero solutions. If X(l)>O, then (1) has the zero solution and its multiplicity is X(1). Next, if the diagram contains two or more points, the non-zero solutions are sought. Having blocks X and Y, we isolate in turn the links of the Newton step-line and find the corresponding pairs (e,R). Assume that, at some step, we have finished considering the next link and have constructed the part of the step-line from a point (X(1), Y(i)) to its vertex To isolate the next link, we find the numbers M,M>K, for which the fraction (X(K), Y(K)). Among these numbers we choose the great[Y(K)-Y(M)]/[X(M)-X(K)] takes its greatest value. est (call it M’); the segment jointing points (X(K),Y(K)) and (X(M’),Y(M’)) is the required link. We find for this link the slope E as an irreducible fraction P/Q. Using blocks X and Y, we discover the points of Newton's diagram which lie on this link. We form the type (12) polynomial and with the aid of a subroutine for evaluating the roots of a polynomial, we find its roots and their multiplicities. This concludes consideration of the link. If the stepline is not yet fully constructed, we proceed to construct the next link. If we are interested in the real solutions of (l), we have to use a subroutine for finding the real roots: if we aim to find all the complex roots of (l), we have to include in the main program a subroutine giving all the complex root. The structure of the main program then remains unchanged regardless of the choice of subroutine. (If the coefficients of the polynomial f(z,W) are rational numbers, and we add a subroutine for calculating the rational roots, we can find only the solutions w(z) for which all the coefficients in expansion (2) are rational numbers.) Before tuning to the next stage of the program. if (1) has an s-tuple solution w(z)=0 f(z, W) by w‘. This amounts to a shift of s cells leftwards of the content (X(1)=5), we divide of block FL (the zeros contained in the first s cells are then lost); block Z and A are not affected. The simplicity with which the division of /(z,w) by &is performed is a further indication of the advantage of our method of writing the polynomial in the computer memory.

Stage 4: change of variables. Every triple (P, Q,R) obtained defines the first term RZ”Q of expansion (2) of the solution of (1). To find the subsequent terms of (2) we have to perform for each triple (P,Q, R) in accordance with Newton's method the change of variables (13). If there is only one such triple, the change (13) can be performed without worrying about sorting the initial polynomial in the memory. If there are several triples, then, before making the change of variables (for one triple), we have to make a copy of the polynomial in the computer memory; then, when computation with the first triple has ended, we can make the change of variables (13) for the other triples. On making the change (13) for the w,(z,);on applying Wewton's first triple, we arrive at a type (1) equation in the function defining the second term in algorithm to the new equations, we obtain the triple (P’,Q’,R’), expansion (2). And again, if this triple is not unique, we have to make a copy of the second polynomial. The process of computation starts to branch, and it can clearly be written in the form of a tree.

polynomial

polynomial Fig.4

Let u be the number of cells employed Copies of the polynomials are made as follows. copy of Then, to obtain a in each of blocks 2 and A in writing the polynomial f(z,w). the content of each cell used in writing it is duplicated in cells with numbers f(G W), increased by U units (in other words, the copy is obtained by shifting lJ cells rightwards the TO COPYblock FL, we use the extra block ELM (Fig.4); after contents of blocks Z and A). obtaining in block FLM the copy of block FL, the content of each cell of block FL is increased This, in the first cell of block FL there appears the number of the pair of cells by U units. of blocks 2 and A, in which is written the last monomial of f(Z,W). containing W to the zeroth The situation is similar as regards the other cells of block FL. degree, In short, the copy of f(z,w) is written in the three blocks Z, A, FL in the same form

10

as that in which the polynomial f(z,m) itself is previously written. of variables (13) is made in the copy. This is done in two stages:

After

this, the change

a) ~-z,~,w-z,% ; only the elements of block 2 are transformed; b) u--w,+R (the idependent variable is not changed); this replacement is made, using the previously prepared matrix of binomial coefficients, calculated by Pascal's triangle method. In short, taking one of triples (P,Q,R), we copy our polynomial and make the change of variables in the copy. We then return to stage 2. We construct the Newton diagram, study it (stage 3) and find the triple (p',Q',R') for the new polynomial. If there is more than one triple, the polynomial obtained after the change of variables (stage 4), is again copied in blocks 2 and A (shifting it rightwards). Thus we now have three polynomials (Fig.4) located in turn in blocks 2 and A. The new block FL is copied in block FLM behind the old block FL already present. We proceed thus until the required number of terms in expansion (21 has been computed. This number is either specified in advance, or else it is found from some accuracy condition. In actual computations the root R is often simple. Then, the triples (p',Q', R'), (P", o", R"),... in expansion (2) are found uniquely, each coefficient R', R”,... being found from a linear equation (see proposition 3). In this case the computations can be simplified, and provision is accordingly made in the program for bypassing certain blocks (e.g., the polynomial copying block, or access to the polynomial root-finding subroutine). Note. Only one-dimensional blocks are used in the program; this assists in speeding the program and economizing on memory. It has been mentioned that the polynomial is more economicallywrittenina one-dimensional rather than two-dimensional block. If we introduce a twodimensional block for copying the polynomials and each new polynomial obtained is written in a new row of this block, this likewise requires a superfluous amount of memory, since the number of monomials in the polynomials obtained during the computation in general increases considerably (see (13) and Fig.4), and in addition, the number of polynomials which simultaneously have to be stored in the memory is not known in advance. The distribution of the memory in the program considered here is in fact dynamic.

4.

On a BESM-6 computer

we solved the equation. (14)

u (0) =o,

5zw+2z’+2s~w+2rd+2~~*-~~-2~w’-o, and obtained

the solutions w=-0.5z+0.1252-0.0625z’+0.03906251’-0.027343~5z’+ o.0205078iz*-0.016113282’+0.01309204z’-0.~19io34z~~ loe-~+~~-&‘+16z’-32z’+64z4-i28z’+256z1-5i2z0~

These solutions

are correct,

since

(14) is equivalent

to the equation

[l+ro-(l+z)“~‘][l+w-(1+2z)-‘]=0. (14) is approximately 1 set, and the program

The time taken to solve translation time is 35 sec. 5. Consider a generalization of the method of Newton diagrams, useful for solving differential equations of any order (including those not solved for the highest derivative). Suppose we are given the equation F(z, w, w',...,d’))=O, where w'~'-w'~'(z),p-0, i,...,8, is the p-th derivative be expressible as

F (z, w, a,. . . , u.) a=&...,

The function

P is assumed to

A,,(Z)W”_‘%’ m-4

where

(15)

of w(z).

a.) is a multi-index,

( 16)

,m,<1

and all coefficients

A,’

and powers

p..' are real

. . . (in particular, F may be an analytic function). Suppose that numbers, A,_*+O, &‘
where O
4

.... II.

w(z)of

(15), written as the series u)(z)=Rz‘+R’za’+ ..,

we can use the following

Step

1.

On the

is replaced by the expression

as a result we obtain

algorithm,

the proof of which

(18) is given below.

right-hand side of (151, every derivative w'p',p=~, 2,..., (e is a parameter to be defined later); a(~-1) . . . (e-p+i)w/z’

the function z,w,ez

z

,...,

e(e-I)...

!19)

11

tr,

(‘4,0z~:=+ . . .) Wm-‘-’(e:)a‘...[e...

(e-s+l):]9=

2

w"[a.p(e)2'~+a,,(e)zT"~+...],

n-0

n-0 i.,
where rflO<~",<. .... &,(e)#O, j=O, I,..., and anj(e) are polynomials (here and below, we assume that n only runs over the integers for which the expression withw"in the brackets of the last sum is not identically zero).

step 2. For the function f(z,lo;e), dependent on the variables z, w and on the parameter we construct Newton's diagram fit will not depend on e, since the power ~~0 does not er depend one ), and we choose a number m in accordance with (S), where it can be assumed without loss of generality that y,,,&O; we plot the points (n, reO),r&m, in the coordinate plane. The absence of a point of the Newton diagram on the ordinate axis implies that the function is a solution of (15). We then draw the Newton step-line L. With every point (n, w(z)=0 T-0) of this step-line we associate a polynomial a.,(e) (Fig.2). Step 3.a. For every link kof the step-line we find its slope; the parameter e is put equal to this slope, while the R, R#O, in expansion (181, corresponding to link L, , is found from the equation (20) a,,,(e)R<+...+a,,.(e)R'r=O, where i,,..) i, are the abscissae terms in expansions (18). b. For each point (n, y,), positive roots of the equation

of points which

lying on Li.

is a vertex

The pairs

of Newton's

(e, R)

will define

step-line,

we evaluate

the first the real 123.)

a..(e)=O.

For each root e, we construct the line with slope e, passing through the point (n,~,~) see Fig.2; if the entire step-line (excluding the points (n, Y.~)) lies strictly above this straight line, then, for any non-zero number R, the pair (e,R) will define the first term of expansion (18). (Here, eis not necessarily equal to a rational number even when F is an analytic function). Thus, the distinct from the algorithm for finding implicit functions, the pairs here, defining the first term in expansions (18), may appear in two ways (see Paras. a and b). Notice that the left-hand sides of (20) and (21) are polynomials, so that the solution of these equations is greatly simplified. Step 4.

For each of the pairs

(e, R), obtained

at step

3, we substitute

into

(15) (22)

tu-z'(R+w,);

w,(z) ; we then return with the result that a type (151 equation is obtained for the function to step 1. As a result of the operation of algorithm II, arbitrary constants may appear in the solutions m(z)(see Step 3. b and step 3.a, when all the coefficients of polynomial (20) are s=l in spite of the presence zero). Notice than an arbitrary constant may also appear when of initial condition (17) (the equation w'z-w=O has the solution w=Cz, -ml (the equation &'z+ w-z). d-1=0, w(O)-0. has the unique solution If, apart from (15) and condition (17), we are also given the initial conditions U'P'(O)-Es,

p==l,2,...,

s-l,

(23)

then instead of seeking the general solution of (15), satisfying (171, and then finding the values of the arbitrary constants, ensuring satisfaction of conditions (23) (which is often not easy to do), we can find directly the solutions of (15), satisfying (17) and (23). For this, the change of required function 1

u== w-bzy-&zZ-...[

(s_l)l Liz’-‘] -$ A_.-

has to be made in (15), and then the result solved for u(z) with the single condition u(O)=O. Similarly, on changing the required function, we can solve problem (15)-(17) under some other auxiliary conditions. We now consider the proof of algorithm II. Let some function m(z), written as series (18), be a solution of (15). We shall show that this series is then among those resulting from the application of algorithm II to (15). Take first the case when series (18) for w(z) is not identically zero, i.e., all the coefficients R, H',... , can be assumed to be non-zero. BY the theorem on the differentiationofa convergent power series, we have W'(Z)-sRz'-'+e'R'z"-'+ ...=eRz'-'+o (z’-‘).

Since R+O. then w(z)-Rz’, Z-CO (i.e., lim[w(z)/(fiz’)l=i), and hence z'-'-wl(Rz) , z+O ). Consequently, w'(z)-ew/z+o(w/z) as z-0. On differentiating series (18) p times, p=l, 2,..., s, we similarly obtain an expression for the p-th derivative of the function u)(z): ~'~'(I)_s(s-i)...(s-P+i)$+O(~), Since w(2) is the solution 2,..., z, into the function

z-co, p==l,2,...,s.

(24)

of (151, then, on substituting w(z) for w and u)(p)(z) for u,, p-l, F(z, 1u,Us,...,u.), given by series (16), identical zero must be

12 obtained, and in particular, the sum of the lowest terms must vanish. By (24), all these lowest terms are the same as the lowest terms of series j(z,m(z); S) , see (191, regardless they will certainly be contained among the following terms: of the value of e>O, Qoe(E)ZTc,

a,,(E)Z’“(RZ’),

. . . , amo(E)Z1r(RZ*)“,

where m is given by condition (5) (cf. (8)). The requirement that the sum of the lowest terms in (25) be identically zero gives us the conditions for finding the possible pairs (S.R); here, unlike the method of finding implicit functions, there are two possibilities. condition (10) holds, then the pair A. If e is such that, for certain i, j, iPj, 092,jGnt, (e,R) is found in the same way as in Para.2: e is the slope of a link of the Newton step-line, and R is calculated from (20). B. If e is not the slope of a link, then the lowest term with this eis unique, let the term be It must vanish identically, which is equivalent to condition (2!.). However, not every positive root of this equation can be the power of the first term in expansion (18); it is also necessary for the power ~,,~+ne of the term (26) to be less than the power of any other term in (251, i.e., that the line with slope E, passing through the point (n, 'jno)of Newton's diagram, lie strictly below all the lines with slope e, passing through the other points This is equivalent to the entire Newton step-line, except for of the diagram, see Fig-l. the point (n, T..~),lying above the line with slope e, passing through the point (n. rmr) (clearly, of the points of the diagram, only the vertex of the step-line, see Fig.2, can have this For the e thus found, given any non-zero R, the pair (S,R) will define the first property). term in expansion (18). are performed, we find In short, we have shown that, when steps l-3of algorithm II the pair (e,R), corresponding to a first term in expansion (18) of solution W(Z) . Next, making the substitution (22) in accordance with step 4, we obtain a new equation of type (15), where S,=e'-e, and UJ,(Z) can one solution of which is w,(z); obviously, ~,(z)=R‘z*+o(z’~), the algorithm, we find, be written as a series of type (18); on again performing steps l-3of in particular, the pair (e,,R’), and from it, the pair (e',R'), giving the second term of expansion (181, etc. The pairs (e, R), (e’,R’), ..., are thus constructed, giving the entire expansion (18) of the function w(Z) (if w(z)*O). Now consider the case when w(z)=O. Then, f(z,0, e)==F(z,O,..,, 0)-O, i.e., there is no in the last sum (19). there is no point on the Newton diagram lying on the term with n-0 ordinate axis, and hence step 2 of algorithm II shows that w(z)-0 is a solution of (15). In short, all the solutions of (15) that can be written as a series (181, will be found by algorithm II. However, differential equations of type (15) may have solutions not expressible as a series (18), though there exist for these solutions and all their derivatives of up to and including order s, type (18) series which are asymptotically convergent to the solutions as To be more precise, if, for Algorithm II may also be used to find these solutions. 2+0. the s-th derivative of a solution w(z)of (15) there existsaseries (181, asymptotically convergent to m(')(z) as z+O, then the function ID(Z)and all its derivatives I@(z), p-l, 2, . . . , s-i, can be written as asymptotic series of type (18), and as a result of applying algorithrl II to (15), we construct a series (18). asymptotically convergent to ID(Z) as z+O. For, to prove algorithm II, it is sufficient to satisfy (241, which in turn follow from We know, representation (18) and the theorem on differentiation of a convergent power series. see /5/, that, if there exists for the derivative of a function an asymptotic series of type (18), and this derivative is continuous, then there is also an asymptotic series of type (18) for the function itself; the series for the derivative being obtained from the series for Hence our above assertion holds. the function by term-by-term differentiation. A differential equation of type (15) may have a solution for which there is no series (18) convergent to it even asymptotically as z-+0. For instance, the equation -&-z‘(uJ')'w(z)=zsin (I/Z-!-~),where C is an 2z'w'w+w'+zzw*=O, w(O)=O, has, in particular, the solution While it is quite impossible to describe these "poor" solutions by means arbitrary constant. of power series, this is hardly of practical interest. However, Computations by algorithm II may be made either "manually," or by computer. there are difficulties in writing a program suitable for solving (15) of general type, due If we confine ourselves to the class to the arbitrary constants that make their appearance. of equations whose solutions do not contain arbitrary constants (this may be known e.g., from physical considerations), or else to the solutions which do not contain arbitrary constants, then no difficulties arise, and the powers e,e'.... for the analytic function F are rational numbers (the case when the solutions do not involve arbitrary constants might be called the of the initial data). case of "pure branching," i.e., branching involving no indeterminacy Then, for a computer realization of algorithm II, it is best ot combine the performance of steps 1 and 2, formula conversions with unknown parameter 6 are then unnecessary, and the program for solving differential equations can be written in a similar way to that considered For writing function F it is best to introduce into above for finding implicit functions. the memory two one-dimensional blocks z and A (in the same way as in Para.3) and a two-dimensional block W, one dimension of which is s, and the other is then length of blocks Z and A, and in addition (just as in Para.3), a small one-dimensional block of powers FL. Comparing Fig.5 with Fig.3, we can see how the information about the left-hand side of (15) is written: in the cOluSinS of Fig.5 with numbers from FL(n)+% to FL(n+l) inclusive is written the sum (see (19))

13

where OC~~N,N is the degree of polynomial F in the set of variables UI, w’....,w(~) (a computer realization of algorithm II works only with finite segments of a power series, i.e., with polynomials); FL(O) is assumed to be zero. I

n n+1 M+!

Fig.5 To be more precise,

if the number

J

satisfies

the condition (27)

FL(n)+f
then in the J-th cells of blocks 2 and A, and in s cells of block W, W(I,p). p=l,Z,....s, is located the information about the next term of the expansion P(r,lo.w',...,u)(')) in series (16): *(,)wt”-r’c’.

After replacing t(z.m;e):

oc.

+wJ.~)I)(~‘)~v.

II ...

ur(p)-e(e-l) ...(e-p+tj~/~p.this term becomes

(Utw)“,Js

qzzv,,

the following

term of the series

*(,)w"lz',,-~l,.',L--. -w,I..,‘~W(,.‘)...[~(E_,)..,(e_l+l)]w(l..,.

(28)

Since we are assuming that no arbitrary constants appear, the pairs (e.R) are simply found __ in accordance with step 3.a. For each n,O=ZnCN, among the numbers J satisfying condition (27). we choose and fix in the memory the J for which * (29) Z(J)W(J, P)P

C. p-1

takes its least value Ino. Knowing 1-0, we can find the number m from (S), construct the Newton diagram, and draw the Newton step-line. Then, for each link L{ of the step-line, we perform the following actions: find the slope e of the link, find the numbers n for which lies on link Lo , and for each such n choose a number J, satisfying condition point (4~4 (27), with which the quantity (29) is equal to I"o (these numbers have already been distinguishedand stored), then for all such J we compute (see (28)) A(J)cw~'~l)...[e(c-i)...(e-a+l)]w"~"; on adding all these expressions, corresponding to fixed n, we find the value of the polynomial Q(8), write (20), and by solving it, obtain all the R corresponding to the slope e. After this, for each of the pairs (e,R) obtained we make the substitution (13). Notice in conclusion that algorithm II is also important because a number of mathematical problems reduces to solving differential equations. For instance, non-elementary functions can often be specified by means of an Eq.(15) and a condition (17). It is sometimes possible to reduce's system of differential equations to a single equation of higher order. Finally, in some cases Eq.(l), where f is an analytic function, can be reduced to an Eq.(15), where F is a polynomial, so that no difficulties arise concerning the expansion of the function f in series (2) and the replacement of this series by a finite segment of it, i.e., a polynomial. For instance, it is easily shown that (w')'(r'-l)+ 1-o z-sinu, * 1 w(O)- 0

co(O)=0 ( m'(O)-i

If we now pass to above modification The author

(15), which is equivalent to (l), no arbitrary constants appear, and the of algorithm II, suitable for computer realization, can be used.

sincerely

thanks N.N. Moiseev

for his guidance.

REFERENCES 1. MOISEEV N.N., Mathematical analiza), Nauka, 1981.

problems

of systems analysis

(Matematicheskie

zadachi

sistemnogo

14 2. KAPASNOSEL'SKII M.A., et al., Approximate solutions of operator equations (Priblizhennye resheniya operatornykh uravnenii), Nauka, Moscow, 1969. 3. VAINBERG M.M. andTRENOGIN V.A., Theory of branching of solutions of non-linear equations (Teoriya vetvleniya reshenii nelineinykh uravenii), Nauka, Moscow, 1969. 4. CHEBOTAREV N-G., Theory of algebraic functions (Teoriya algebraicheskikh funksii), Gostekhteorizdat, Moscow-Leningrad, 1948. 5. OLVER F., Introduction to astmptotic methods and special functions (Vvedenie v asimptoticheskie metody i spetsial'nye funksii), Nauka, Moscow, 1978. Translated

U.S.S.R. Comput.Maths.Math.Phys.. Printed in Great Britain

by D.E.B.

0041-5553/84 $lo.cO+o.Oo 01985 Pergamon Press Ltd.

Vo1.24,No.4,pp.l4-18,1984

TWO SCHEMES FOR A NON-LINEAR METHOD OF OPTIMIZATION IN EXTREMAL PROBLEMS* A.A. TRET'YAKOV Two schemes of the method of non-linear descent are proposed for solving the unconstrained minimization problem q(r)*min,SPE,,Z the method with independent choice of coefficients, and a method of quadratic type. Conditions are stated for the existence of a minimizing arc for the general problem of mathematical programming. The methods considered below for solving both unconstrained and constrained extremal problems have a scheme which differs from the traditional linear scheme of iterative optimlzation methods. Most methods are constructed using the recurrence relation (1) tr+,=xr+aksr, where sris the direction of descent, chosen by a certain rule, c& is the length of the step in this direction, and zk is the next approximation. The scheme in the present paper is (2) Z~+,-k+~d~+~&, where sh,d, are directions of descent, or the minimizing pair, which specifies the arc of descent, and a*, b& are the step lengths in these directions. The minimization in this scheme is not along a straight line, but along a certain set, specified by directions s*,d, and coefficients In this sense the process is non-linear with respect to the minimization oL*, Bk. set at each step. Construction (2) is considered for the following reasons: to solve problems with strong ravines, where linear methods usually converge slowly or not at all, due to the influence of rounding errors, to solve problems of degenerate type, when the first derivatives of the objective function or the functions of constraints vanish at the minimum point, and to solve problems of non-linear programming with an admissible set of complex configuration. Below we propose two versions. of the methods of scheme (2) and prove convergence theorems. 1.

Non-linear methods of minimization with independent choice of coeffi-

cients. The unconstrained

extremal

problem cp(z)+min, ZEE, will be solved by a method whose scheme is given by the relations (3a)

x~+,-z~-~r-_B~& (cp'(Sh), SJW i.e., & and&are directions from the condition

$'(zr),&)>O,

of decrease. q (z,-aa8k-Adk)

(3b)

Il~~ll=lldrll==~,

It would be natural

to seek coefficients

&and

p*

= nf,n q b-~,-BdJ.

In general, however, this problem can only be solved approximately, which causes further difficulties. Henceforth it will always be assumed that the function v(z) belongs to the class C','(X,), X,=(r~E.l~(z)scp(zo)), and we shall seek coefficients QI, p,, satisfying the condition 'P(.~~)~(z~+,)~Q~Q~(~'(z;), rr)+qrfUcp'(zl), d&-La&

(41

where L is the Lipschitz constant of the gradient of v(z), and ql is a coefficient. prove below the existence of a,,p,, guaranteeing that (4) holds. We have: cients

Lemma. ah, fl,

Let a sequence {z,) be constructed satisfy condition (4). Then,

accordingtoscheme

'~(fr)--(~(~~+,)>~-'q~Ar'll~'(~~)II*, where *Zh.vychisl.Mat.mat.Fiz.,24,7,986-992,1984

(3) and let the

We shall coeffi(5)