Solution on a digital computer of the finite operator equation Ax=? by the method of continuous steepest descent
Abstract Generated abstract
This paper proposes a digital-computer method for solving finite, generally nonlinear operator equations of the form Ax = 0. It reformulates the problem through the residual norm functional and derives a system of ordinary differential equations whose trajectories follow the normals to nested level surfaces, interpreted as paths of continuous steepest descent toward a root. The method is presented as a simpler alternative to broken-line descent schemes and as analogous to analog-computer approaches, but intended for higher accuracy on digital machines. Because the differential system becomes indeterminate near a root, the paper recommends using a relatively large-step numerical integration procedure for the descent phase and switching to Newton’s method once the residual functional is sufficiently small.
Full Text
SOLUTION ON A DIGITAL COMPUTER OF THE FINITE OPERATOR EQUATION \(Ax=\Theta\) BY THE METHOD OF CONTINUOUS STEEPEST DESCENT
B. F. Krutoi
(Presented by the scientific-technical seminar of the Department of Applied Mathematics and the Educational Computing Laboratory)
Solution of the equation
\[ Ax=\Theta \tag{1} \]
with a given, generally nonlinear, finite operator \(A\), the zero vector \(\Theta=(0^1,0^2,\ldots,0^n)\), and the sought vector \(x=(x^1,x^2,\ldots,x^n)\), is usually carried out on a digital computer by broken-line descent methods [1], whose algorithms are rather complicated and do not always possess sufficiently good convergence. Meanwhile, problems of this kind have long been solved approximately on analog computers by reducing them to ordinary differential equations with independent time \(t\) [2]. Therefore, in the present article, for the solution of the same finite equation (1), a simpler and more reliable method will be considered: not broken-line, but continuous steepest descent. The proposed steepest-descent method, carried out on a digital computer, is somewhat reminiscent of the known methods for solving equation (1) on an analog computer, but differs from them in its higher accuracy and requires, near the root \(x\), a transition to Newton’s method. The essence of the method of continuous steepest descent is as follows.
Substitute into equation (1), in place of the sought \(x=(x^1,x^2,\ldots,x^n)\), an arbitrary vector \(\xi=(\xi^1,\xi^2,\ldots,\xi^n)\). Then equation (1) takes the form
\[ A\xi=r, \tag{2} \]
where \(r=(r^1,r^2,\ldots,r^n)\) is the residual vector. We then form the positive-definite functional \(\Phi(\xi)\)
\[ \Phi(\xi)=(A\xi,A\xi)=(r,r)\ge 0, \tag{3} \]
which becomes 0 when \(\xi=x\). If in this relation the arbitrary vector \(\xi=(\xi^1,\xi^2,\ldots,\xi^n)\) is replaced by the given vector \(\xi_0=(\xi_0^1,\xi_0^2,\ldots,\xi_0^n)\), then the equality
\[ \Phi=\Phi(\xi)=\Phi(\xi_0)=\Phi_0=(r_0,r_0), \tag{4} \]
where \(r=(r_0^1,r_0^2,\ldots,r_0^n)\), gives the equation of a level surface passing through the point \(M_0=M(\xi_0)\) of the geometric space \(D^{(n)}\).
Let us now draw, through an arbitrary point \(M=M(\xi)\) of the level surface \(\Phi=\Phi_0\), the tangent plane and the surface normal. Let \(d\bar r=d\xi^i \bar e_i\) be an arbitrary infinitesimal displacement from the indicated point \(M\), lying in the tangent plane. Then the condition is valid
\[ \bar N \odot d\bar r=0, \tag{5} \]
expressing the mutual perpendicularity of the surface normal \(\bar N\) and the tangential displacement \(d\bar r\). If we specify the surface normal \(\bar N\) by the expansion \(\bar N=N_\lambda \bar e^\lambda\), then condition (5) takes the form
\[ \bar N \odot d\bar r = N_\lambda \bar e^\lambda \odot d\xi^i \bar e_i = N_\lambda d\xi^i(\bar e_\lambda \odot \bar e^i) = N_i d\xi^i, \tag{5*} \]
where we have taken into account the fundamental relation in tensor analysis
\[ \bar e_\lambda \odot \bar e^i=\delta_\lambda^i= \begin{cases} +1, & \text{if } \lambda=i,\\ 0, & \text{if } \lambda\ne i. \end{cases} \tag{6} \]
On the other hand, since the left-hand side of equality (4) \(\Phi(\xi)=\Phi(\xi_0)\) does not change when different numbers \(\Phi(\xi_0)\) are substituted into the right-hand side, from both sides of this equality one may take the total differential \(d\), and we obtain
\[ d\Phi=\frac{\partial \Phi}{\partial \xi^s}d\xi^s=d\Phi_0=0. \tag{7} \]
Comparison of (5*) and (7) shows that
\[ N_s=\frac{\partial \Phi}{\partial \xi^s}. \tag{8} \]
Hence follows the following refined form of condition (5):
\[ \bar N \odot d\bar r = \frac{\partial \Phi}{\partial \xi^s}d\xi^s = 0. \tag{5**} \]
Next we shall construct the vector differential equation of the normal \(\bar N\) to the level surface \(\Phi(\xi)=\Phi(\xi_0)\) at an arbitrary point \(M(\xi)\) of this surface. If by \(d\bar r\) we now denote an infinitesimal displacement from the indicated point \(M(\xi)\), directed along the normal \(\bar N\), then we obtain the following vector equation of this normal:
\[ \bar N = N_s \bar e^s = \frac{\partial \Phi}{\partial \xi^s}\bar e^s = \mu\, d\bar r = \mu\, d\xi^i \bar e_i, \tag{9} \]
where \(\mu\) is an arbitrary numerical factor. To pass from here to the corresponding differential equations in the components \(N_s\), \(d\xi^i\), and to obtain the simplest form of such equations, suppose that the vector \(\xi\) is given in Cartesian rectangular coordinates. Then
\[ \text{1) } \bar e_i=\bar e_i=\bar e^i=\bar e^i,\qquad \text{2) } d\bar r=d\xi^i\bar e_i=d\xi_i\bar e^i,\qquad \text{3) } d\xi^i=d\xi_i, \tag{10} \]
\[ 4)\quad w_{\alpha\beta}=\bar e_\alpha\otimes \bar e_\beta=\bar e^{\alpha}\otimes \bar e^{-\beta}=w^{\alpha\beta} = \begin{cases} w_{\alpha\alpha}=w^{\alpha\alpha}=+1, & \text{if } \beta=\alpha,\\ 0, & \text{if } \beta\ne\alpha. \end{cases} \]
Hence, taking into account the relation \(N_s=\partial\Phi/\partial\xi^s\), we find the following representation of vector equation (9) in the components \(N_s,d\xi^i\):
\[ 1)\frac{\partial\Phi}{\partial\xi^s}=\varphi_s(\xi),\quad 2)\quad \frac{\varphi_1(\xi)}{d\xi^1}= \frac{\varphi_2(\xi)}{d\xi^2}=\cdots= \frac{\varphi_n(\xi)}{d\xi^n}=\mu. \tag{11} \]
In the chain of \(n\) links (10.2), let us take \(\xi^1\) as the independent variable, and regard the remaining \(n-1\) variables \(\xi^s\) as dependent. Then we arrive at the following system of \(n-1\) first-order differential equations:
\[ \frac{d\xi^s}{d\xi^1}= \frac{\varphi_s(\xi)}{\varphi_1(\xi)} =\psi_s(\xi),\quad (s=2,3,\ldots,n), \tag{12} \]
which determine the direction of the normal \(\bar N\) at an arbitrary point \(M(\xi)\) of the level surface \(\Phi(\xi)=\Phi(\xi_0)\).
Suppose now that in the geometric space \(D^{(n)}\) we have constructed a continuous set of nested level surfaces \(\Phi_i(\xi)=\Phi(\xi_i)\), \((i=0,1,2,\ldots)\), passing through selected points \(M_i=M(\xi_i)\) of the space \(D^{(n)}\) and satisfying the condition \(\Phi_i(\xi)>\Phi_{i+1}(\xi)\). In this case the Cauchy problem
\[ 1)\frac{\partial\Phi(\xi)}{\partial\xi^s}=\varphi_s(\xi),\quad 2)\frac{d\xi^s}{d\xi^1}= \frac{\varphi_s(\xi)}{\varphi_1(\xi)} =\psi_s(\xi), \tag{13} \]
\[ 3)\quad (\xi^s)_{\xi^1=\xi_0^1}=\xi_0^s \]
defines a certain continuous thread \(H^{(n)}_0\) in \(D^{(n)}\), which passes through the chosen initial point \(M_0=M(\xi_0)\) of the space \(D^{(n)}\) and pierces transversely each surface \(\Phi_i(\xi)=\Phi(\xi_i)\) from the above-mentioned continuous family of them, subordinated to the condition \(\Phi_i(\xi)<\Phi_{i+1}(\xi)\). A continuous thread \(H^{(n)}\) in \(D^{(n)}\) satisfying such conditions will be called a thread of the shortest continuous descent. We may even assume that the points \(M_i=M(\xi_i)\) of the nested level surfaces \(\Phi_i(\xi)=\Phi(\xi_i)\) are located precisely on the sought thread of shortest continuous descent. We shall use this remark below.
From the chain of reasoning we have given, it now becomes, as it were, clear that the search for roots \(x\) of equation (1)
\[ Ax=\Theta \]
can be reduced to solving the simpler Cauchy problem (when performed on a digital computer) for various initial points \(M_0=M(\xi_0)\). Having chosen one such initial point \(M_0\), we construct, from the intersection points \(M_i\), the thread \(H^{(n)}_0\) of shortest continuous descent issuing from \(M_0\), each time computing the corresponding value \(\Phi(\xi_i)\). That point \(M=M_k=M(\xi_k)\) for which \(\Phi(\xi_k)\approx0\) will be one of the roots \(x\) of equation (1).
In the same case when, instead of a point \(M(\xi_k)\) with mark \(\Phi(\xi_k)\approx0\), we fall into a pit with bottom mark \(\Phi(\xi_k)>0\), it is necessary to choose as ...
initial point \(M_0=M(\xi_0)\) some other point, from which we begin a new continuous descent along the steepest path.
We solve the Cauchy problem (13), of course, numerically, with some step \(h=d\xi^1\). However, in this way we shall not obtain exact values \(\xi_k\) for the roots \(x\) of equation (1) in a finite number of steps. This is due to the fact that, for the root \(x=\xi_k\), the functional \(\Phi(\xi)_{\xi=x}=0\), and therefore in (13)
\[ 1)\ \frac{\partial \Phi(\xi_k)}{\partial \xi^s}=\varphi_s(\xi_k)=0,\qquad 2)\ \psi_s(\xi_k)=\frac{\varphi_s(\xi_k)}{\varphi_1(\xi_k)}=\frac{0}{0}. \tag{14} \]
Consequently, near the root \(\xi_k=x\) the right-hand sides \(\psi_s(\xi_k)\) of equations (13.2) become indeterminate. Therefore, as soon as for some \(\xi_i\) the functional \(\Phi(\xi_i)\) becomes sufficiently small, for example \(\Phi(\xi_i)=0.5\text{–}0.7\), one must switch to Newton’s method. It follows from this that the numerical solution of the Cauchy problem (13) can be carried out with a large step \(h=d\xi^1\), using Euler’s method or a simplified Runge–Kutta method of second-order accuracy.
As for the particular form of the equation \(Ax=\Theta\) being solved, it may be either a single equation \(f(z)=0\), or a system of such equations
\[ f_i(z_1,z_2,\ldots,z_n)=0,\quad (i=1,2,\ldots,n). \]
REFERENCES
- E. D. Booth. Numerical Methods. Moscow, GIFML, 1959.
- V. V. Smolov et al. Continuous-Action Computing Machines. Moscow, “Vysshaya shkola,” 1964.