An Approximate Solver for Symbolic Equations

1 Introduction

Research in symbolic equation solving has focused on ex­ act solution methods. The resulting programs, such as MACSYMA (Mathlab Group, 1983) and PRESS (Bundy and Welham, 1981], either return an exact solution or fail. Yet , scientists and engineers routinely must solve prob­ lems that have no exact closed-form solution. They need an equation solver that finds adequate approximate solu­ tions to such problems, rather than failing. In fact, they often prefer a simple approximate solution even when a complicated exact solution is available. For example, the equation x 5 + x = k cannot be solved for x in closed form, but x = k is an accurate approximation to the solution for k near 0. One might well prefer this approx­ imation to the full solution of x 4 + x = k, which is very long and complicated. This paper describes an equation solver, called NEWTON, that finds approximate solutions to parameterized equations in one unknown.

Numeric analysis provides many algorithms for solv­ ing individual equations approximately, but says little about parameterized equations. One could solve numer­ ically for specific parameter values and interpolate the results. This would provide little general understanding and require prohibitive amounts of computation, espe­ cially for equations containing multiple parameters. In­ stead, NEWTON constructs a single parameterized solu­ tion, such as x = k above, for all legal parameter values. Users can instantiate a solution with multiple parame­ ter values, rather than rederiving it numerically for each value. They can examine the influence of parameters on a solution analytically instead of by experimentation.

For example, the approximate solution x = k increases linearly in k.

2 The Algorithm

NEWTON takes as input a parameterized function f(x), an interval domain for x, and a set of constraints. Each constraint is a strict or nonstrict inequality between real­ valued functions of the parameters, for example k > 1 and k 2 :::; 2m. The constraints and /( x) must be ex­ tended elementary functions: polynomials and compo­ sitions of exponentials, logarithms, trigonometric func­ tions, inverse trigonometric functions, and absolute val­ ues. NEWTON assumes that f(x) is differentiable and has a single, simple root 1 in the domain of x. I discuss multiple and nonsimple roots in the concluding section.

NEWTON narrows in on the root in stages. First, it brackets the root within an interval on which f(x) is monotonic by the following algorithm:

1. let [a, b] be the domain of x
2. iff' (x) has a fixed sign on [a, b] then return [a, b]
3. let m =(a+ b)/2; if f(a)f(m) :::; 0 then b +- m else a+- m
4. go to step 2.

If a = -oo or b = oo, it chooses a very small or very large finite number instead. It uses the BOUNDER inequality prover [Sacks, 1987a, Sacks, 1987b] to test whether the inequalities in steps 2 and 3 hold for all parameter val­ ues that satisfy the input constraints. The simple root assumption ensures that the algorithm terminates. For example, the function f(x) = x 4 +kx-16 with parameter k and constraints -1 < k < 1 is monotone on the domain x E [1,10) because f'(x) = 4x 3 +k 2:4-1 = 3 > 0.

Next, NEWTON calculates an initial approximate root, x0 , for f(x). If f(x) is a sum, it finds the addend of max­ imal absolute value among those containing x, deletes the remaining addends that contain x, and solves the resulting equation. 2 In our example, the addends x 4 and kx contain x with lx 4 1 > lkxl, so NEWTON solves x 4 - 16 = 0 to obtain x0 = 2. It rejects x0 = -2 be­ cause that root lies outside the domain. NEWTON tries to calculate x0 for general equations by reducing them to sums, as shown in Table 1. These reductions suf­ fice for all rational functions and for all the examples in this paper. If NEWTON cannot reduce a function or cannot solve a reduced equation, it sets x0 to be the middle of the monotone domain of x. One such case is f(x) = x"'- 2.

1 A root r is simple if f(r) = 0 and f'(r)-:/:. 0.

2 McAllester [McAllester, 1981] uses a similar strategy.

NEWTON tries to improve upon x0 by applying New­ ton's method symbolically. It tests whether the first Newton iterate

is closer than x0 to the root, r, of f(x), that is whether

It must test this condition indirectly because r is un­ known. If the test succeeds, NEWTON assumes that Xi approximates r sufficiently well and returns Xi· This assumption normally works well. In our example,

has a relative error of at most 0.0015 for all k that satisfy the input constraint. A more sophisticated implementa­ tion could check xi against an accurate numeric solution for several parameter values and calculate further iter­ ates when xi proved inadequate. The stopping condition would involve a tradeoff between the accuracy and the complexity of the approximate solution, both of which grow with each iteration.

To understand how NEWTON tests xi, suppose first that f(x) increases monotonically and that x0 > r. Fig­ ure 1 depicts the geometric interpretation of Newton's method in this situation: xi is the intersection point between the tangent to f(x) at x0 and the x axis. As f'(xo) decreases from oo to 0, Xi decreases from x0 to -oo via r. Any xi 2: r satisfies equation (1), as Fig­ ure 1a illustrates. If Xi < r (Figure 1b) , equation (1) becomes

By the mean value theorem, f(x0 ) = f'(~)(x0 - r) for some ~ in (r, x0 ). Substituting into equation (2) and rearranging yields/'(~)< 2f'(x0 ). A sufficient condition for this inequality is that f'(x) < 2f'(x0 ) for all x in [a, x0] . This condition holds in our example:

Figure 2: A sufficient condition for ixo- rl > ixi - ri on a monotone interval [a, b].

for x in [1, 2]. The cases of x0 < r and of f(x) mono­ tone decreasing are analogous. The complete test for equation (1) appears in Figure 2. NEWTON resolves the inequalities with BOUNDER as before.

NEWTON tests progressively narrower intervals until it succeeds or the interval shrinks beneath a prespecified width. It starts with the monotone interval [a, b] and repeatedly bisects the distance from x0 to each current boundary. If the initial guess really is a good candidate for iteration, the shrinking should eventually make the test work, at which point NEWTON returns the first iter­ ate. If not, NEWTON must choose a different value for x0 by bisection and start over. If the test had failed in the example, the next test interval would have been (1.5, 6], splitting the distance from the ends of [1, 10] to x0 = 2.

3 Examples

As another example, consider f(x) = e"' + x- k with domain x E [0, oo) and constraint k ~ 1. The function increases monotonically for all x because f' ( x) = e"' + 1 > 0. It has a unique root because f(O) = 1 - k :::; 0 and f( k) = ek > 0. Approximation is inevitable because the root cannot be expressed in closed form. NEWTON determines that le"' I > lx I and solves e"' - k = 0 to obtain xo = log k . It verifies the convergence condition (Figure 2) and returns

Table 2: Relative error between the root of e"' + x- k and x1 = k log kj(k + 1).

Given the domain x E ( -oo, -1] and no constraints on k, NEWTON determines that lxl > le"'l, obtains x0 = k, verifies convergence, and returns

For f ( x) to have a root, k must satisfy

The relative·error of x1 is 0.2 for this value, 0.007 for k = -1, and less than 0.001 fork< -2.

The dominance argument only requires the sign of k: k > 0 implies le"'l > lxl, whereas k < 0 implies lxl > ie"'l· NEWTON needs stronger constraints, such as those above, to derive dominance. Given only the sign of k, it must calculate x0 by bisection. It can, however, prove convergence and obtain x1 as before.

The previous examples are all separable, that is ex­ pressible in the form f(x, k) = g(x) + h(k). The final example,

with x E (0, oo) and k ~ 0.5, demonstrates that NEW­ TON handles nonseparable functions as well. The func­ tion increases monotonically and has a unique root. The convergence condition holds for all k. Given the modest constraint k < e"'~ - x, NEWTON can calculate x0 = k and

The relative error of x1 is 0 for k = 0.5, 0.048 for k = 1, and less than 0.002 for k > 4. One can obtain even bet­ ter results with a more sophisticated choice of x0 , setting x = 0 in log(x+k), the dominated term of f(x), and solv­ ing the resulting equation to obtain x~ = Jk 2 - log k. NEWTON does not implement this strategy. Newton's method appears to work for 0 < k < 0.5, even though the convergence condition fails. A more experimental implementation could validate x1 empirically on a large number of sample parameter values, rather than give up.

Table 3 summarizes the examples in this paper.

4 Conclusions

This paper describes NEWTON, an equation-solver that finds approximate roots to parameterized functions by a symbolic version of Newton's method. If this fails, it turns to a symbolic bisection algorithm. NEWTON re­ quires that the function have a single, simple root in its domain. One can .handle a nonsimple root, r, of f(x) ·by applying NEWTON to ~t(x) = f(x)ff'(x), which has a simple root at r. One can find multiple roots by parti­ tioning the domain into intervals on which f(x) is mono­ tone, hence has at most one root, and applying NEWTON to .each subinterval [a, b] for which f(a)f(b) :::; 0. I de­ scribe a program that calculates the monotone intervals elsewhere [Sacks, 1985].

NEWTON can only find a root r of f(x) for which the sign of f' ( r) is fixed for all legal parameter values. Oth­ erwise, it can neither bracket the root within a monotone domain nor prove that Newton's method converges. For example, the function f(x) = (ax- 1)e"' has a unique root of 1/a for a =F 0, but /'(1/a) = ae 1 fa can be positive or negative. A possible solution is to partition the space of legal parameter values into regions on which J'(x) has a fixed sign and solve for each subset separately. In our example, one can solve for a > 0 and for a < 0. Partition algorithms exist for polynomials [Arnon et al., 1984]; partitioning other functions is a topic for future research.

The strategy of applying numeric techniques symbol­ ically to parameterized problems has many other appli­ cations. Newton's method generalizes to vector equa­ tions. Additional numeric techniques, such as the se­ cant method and fixed-point iteration [Press et al., 1987], also extend directly to parameterized equations. Other domains for symbolic application of numeric algorithms include Markov theory (which I have worked on with

1. Doyle [Doyle and Sacks, 1989]), integration; and dy­ namic systems theory.

Acknowledgem.ents

Harry Mynick of the Princeton Plasma Physics Labo­ ratory motivated me to develop NEWTON. Forman Ac­ ton, Woody Bledsoe, Joel Friedman, Dick Lipton, Harry Mynick, and especially Mike Wellman provided valuable comments on this paper.