





























                              Research Report AI-1989-02

                           Artificial Intelligence Programs

                              The University of Georgia

                                Athens, Georgia 30602


                                         Available by ftp from

                                           aisun1.ai.uga.edu

                                            (128.192.12.9)


                                            Series editor: 

                                           Michael Covington

                                      mcovingt@aisun1.ai.uga.edu










                             A Numerical Equation Solver in Prolog


                                 Michael A. Covington


                                   Artificial Intelligence Programs

                                       The University of Georgia

                                         Athens, Georgia 30602


                                              March 1989



                        Abstract: The Prolog inference engine can be extended

                        to solve for unknowns in arithmetic equations such as

                        X-1=1/X or X=cos(X), whether or not the equations have

                        analytic solutions. This is done by standard numerical

                        methods, but two features of Prolog make the

                        implementation easy: the ability to treat expressions

                        as data and the ability of the program to extend

                        itself at run time.




                  1. The problem


                        The Prolog inference engine can solve for any unknown in

                  symbolic queries, but not in arithmetic queries. For example,

                  given the fact


                              father(michael,sharon).


                  one can ask


                              ?- father(michael,X).


                  and get the answer sharon, or ask


                              ?- father(X,sharon).


                  and get the answer michael. This interchangeability of unknowns

                  extends to complex symbolic manipulations (e.g., append can be

                  used to split a list as well as to concatenate lists), but not to

                  arithmetic. 


                        Prolog handles arithmetic the way Fortran did thirty years

                  ago: the unknown can only be a single variable on the left of the

                  operator is, and everything on the right must be known. Thus


                              ?- X is 2 + 2.


                  gets the answer 4, but


                              ?- X is 1 + 1 / X.









                                                                                   2


                  fails or raises an error condition.


                        The excuse for this restriction is that Prolog cannot search

                  the set of real numbers the way it searches the symbols in a

                  knowledge base. As far as exhaustive search goes, this is true.

                  However, mathematicians have been using heuristic searches to

                  solve equations since the days of Isaac Newton. The procedures

                  given here implement one such method, making it possible to have

                  dialogues with the computer such as:


                              ?- solve( X = 1 + 1 / X ).

                              X = 1.618034


                              ?- solve( X = cos(X) ). 

                              X = 0.739085


                  and so on.



                  2. The solution


                        Prolog is an ideal language for solving equations for two

                  reasons: equations can be treated as data, and the program can

                  modify itself. A procedure can accept expressions as parameters,

                  then evaluate them or even create procedures to evaluate them. In

                  Pascal or C, by contrast, there is no simple way to introduce a

                  wholly new equation into the program at run time.


                        The program given here solves equations by the secant

                  method, which is one of the simplest numerical methods, though not

                  the most robust. A different method can easily be substituted once

                  the framework of the program is in place. Standard (Edinburgh-

                  compatible) Prolog is required; Turbo Prolog programs cannot

                  modify themselves in the necessary way.


                        To solve the equation 


                              Left = Right 


                  the secant method uses the function 


                              Dif(x) = Left - Right 


                  where Left and Right are expressions that contain x. The problem

                  is then to search for an x such that Dif(x) = 0.


                        The search is begun by taking two guesses at x and comparing

                  the values of Dif(x) for each. One of them will normally be closer

                  to zero than the other. From this information the computer can

                  tell whether to move toward higher or lower guesses. In fact, by

                  assuming that Dif(x) increases or decreases linearly with x, the

                  computer can estimate how far to move. (This is why it's called









                                                                                   3


                  the secant method -- given two guesses, the third guess is formed

                  by extending a secant line on the graph of the function.) 


                        Success is not guaranteed -- the two Dif values could be

                  equal, or the estimate of how far to move could be misleading --

                  but the procedure usually converges on a solution in just a few

                  iterations. Listing 1 shows the algorithm in pseudocode form.



                  3. Finding the unknown


                        Expressing this in Prolog boils down to two things: setting

                  up the problem, and then performing the computation. The setting

                  up is done by solve, which calls free_in, define_dif, and

                  solve_for (Listing 2).


                        The first step is to identify the unknown -- that is, to

                  pick out the free variable in the equation to be solved. This is

                  done by procedure free_in, which finds the free (uninstantiated)

                  variables in any Prolog term. This is more general than what we

                  need, but it's always useful to build a general-purpose tool.


                        If the term contains a free variable, there are two

                  possibilities: either the term is the variable, in which case the

                  search is over, or else the term has a variable somewhere in its

                  argument structure. free_in has a clause for each of these cases.


                        The second case requires that the term be decomposed into a

                  list. The built-in predicate "univ" (=..) does this. For example,


                              a(b,c(d),e) =.. [a,b,c(d),e].


                              2+3+X =.. ['+',2,3+X]


                  Even lists can be split this way, because any list is really a

                  two-argument structure with the dot ('.') as its principal

                  functor. That is, the list [a,b,c] is really a.(b.(c.[])), though

                  not all Prologs allow you to write it that way. So,


                              [a,b,c] =.. ['.',a,[b,c]].


                  Thus a list is not a special case; it can be treated just like any

                  other complex term.


                        In free_in, the statement


                              Term =.. [_,Arg|Args].


                  discards the functor, which can't be a variable anyway, and

                  obtains two things: the first argument, Arg, and the list of

                  subsequent arguments, Args. It is then straightforward to search

                  for variables in both Arg and Args. Further, because Arg can be a









                                                                                   4


                  single variable, the first clause has a chance to terminate the

                  recursion; and if it isn't, whatever is the first element of Args

                  on this pass will be Arg on the next recursive pass and will get

                  examined then.


                        There is, however, a special case to rule out. The term [[]]

                  decomposes into ['.',[],[]], givng Arg = [] and Args = [[]], which

                  would lead to an endless loop. For this reason, free_in explicitly

                  tests for this term and rejects it.



                  4. Defining a procedure


                        The next step is to define a procedure to compute the Dif

                  function. Recall that the argument of solve is an equation in the

                  form Left=Right. Clearly, the Dif function is obtained by

                  evaluating Left-Right. But how is this done?


                        There are two possibilities. One possibility would be to

                  pass along the expression Left-Right and evaluate it whenever

                  needed. This is easily done, because


                              X is Y.


                  will evaluate whatever expression Y is instantiated to. 


                        But the other, faster, possibility is to define a procedure

                  to do the evaluating. That's what define_dif does; it creates a

                  procedure such as


                              dif(X,Dif) :- Dif is X - cos(X).


                  using whatever expressions the user originally supplied. The

                  result is a procedure called dif that accepts a value of X and

                  returns the corresponding value of the Dif function. In ALS

                  Prolog, this dif procedure is compiled into threaded code when

                  assert places it into the program; it runs just as fast as if it

                  had been supplied by the original programmer. Other Prologs run it

                  interpretively. 


                        What connects the variable X to the expression Left-Right in

                  which it is supposed to occur? This question addresses the heart

                  of Prolog's variable scoping system. It isn't enough simply that

                  it is called X; like-named variables in Prolog are not the same

                  unless they occur in the same rule, fact, or goal. 


                        That's why so many goals in this program have both X and

                  Left=Right (or Left-Right) as arguments. Initially, free_in takes

                  Left and Right and finds a variable in them. This variable may

                  have any name, but it is unified with the variable X in solve.

                  This same X is then passed, along with Left and Right, to

                  define_dif, which uses it in creating the dif procedure.









                                                                                   5


                  Thereafter, the X in dif is guaranteed to be a variable that

                  occurs in Left-Right, also in dif, regardless of what the user

                  originally called it.



                  5. Solving the equation


                        The last step is to implement the secant method (Listing 3).

                  The pseudocode in Listing 1 undergoes several changes when

                  translated into Prolog. 


                        First, Prolog has no loop constructs, so recursion is used

                  instead. The loop is replaced by the procedure solve_aux, which

                  calls itself. Because the recursive call is the last step of a

                  procedure with no untried alternatives, the compiler converts it

                  back into a loop in machine language, but conceptually, the

                  programmer thinks in terms of recursion.


                        Second, in Prolog there is no way to change the value of an

                  instantiated variable. This means, for example, that there is no

                  Prolog counterpart of


                              Guess1 := Guess2


                  when Guess1 already has a value. Instead, the proper Prolog

                  technique is to pass a new value in the same argument position on

                  the next recursive call. Thus the procedure that begins with


                              solve_aux(...Guess1,Dif1,Guess2...) :- ...


                  ends with the recursive call


                              ... solve_aux(...Guess2,Dif2,Guess3...).


                        Third, there are minor rearrangements to avoid computing

                  Dif(X) more than once with the same value of X. These include the

                  variable Dif2 and the passing of Dif1 as an argument from the

                  previous recursive pass.



                  6. Limits and possibilities


                        This program is intended as a demonstration of the

                  integration of numerical methods into Prolog, not as a

                  demonstration of numerical methods per se.


                        The secant method is simple, but far from perfect. It has

                  trouble with some equations. For example, if two successive Dif

                  values happen to be exactly the same distance from zero, then

                  solve_aux will try to divide by zero. This simply fails in ALS

                  Prolog but may cause an error message in other Prologs. This

                  problem shows up with the equation









                                                                                   6


                              X^2 - 3*X = 0


                  which has Dif=-2 for both of the first two guesses (1 and 2). With

                  a case very close to this, such as X^2 - 3.01*X = 0, we find that

                  although the method should work in principle, in practice the next

                  guess is a long way from the correct solution, and the guesses run

                  wildly out of the range of representable numbers. And with some

                  equations, the guesses will bounce back and forth between two

                  values, not getting better with successive iterations [1].


                        More robust numerical methods can easily be substituted into

                  the same program framework. The ability to solve for more than one

                  unknown is desirable; this could be treated as a multi-variable

                  minimization problem where the goal is to minimize abs(Dif(X))

                  [2]. It is possible to solve systems of nonlinear equations by

                  reducing them to systems of linear equations, which can then be

                  solved by conventional methods.


                        The program was written in ALS Prolog and has been tested in

                  Quintus Prolog. However, other Prologs may require minor

                  modifications. For example, Arity Prolog 4.0 requires spaces

                  before certain opening parentheses (e.g., 2 + (3+4) + 5 rather

                  than 2+(3+4)+5). And it is a general limitation of real-number

                  arithmetic that a negative number cannot be raised to a non-

                  integer power (i.e., 4^2.5 is all right but (-4)^2.5 is not). Some

                  Prologs assume all exponents are non-integer.



                  References


                  [1] Hamming, Richard W., Introduction to Applied Numerical

                  Analysis (New York: McGraw-Hill, 1971), especially pp. 33-55.


                  [2] Press, William H.; Flannery, Brian P.; Teukolsky, Saul A.; and

                  Vetterling, William T., Numerical Recipes: The Art of Scientific

                  Computing (Cambridge: Cambridge University Press, 1986),

                  especially pp. 240-334.









                                                                                   7


                  Listing 1. The secant method algorithm in pseudocode form.



                  To solve Left = Right:


                    function Dif(X) = Left - Right

                      where X occurs in Left and/or Right;


                    procedure Solve;

                    begin

                      Guess1 := 1;

                      Guess2 := 2;

                      repeat

                        Slope := (Dif(Guess2)-Dif(Guess1))/(Guess2-Guess1);

                        Guess1 := Guess2;

                        Guess2 := Guess2 - (Dif(Guess2)/Slope)

                      until Guess2 is sufficiently close to Guess1;

                      result is Guess2

                    end.









                                                                                   8


                  Listing 2. Procedures to set up the problem.



                  % solve(Left=Right)

                          %

                          %  On entry, Left=Right is an arithmetic

                          %  equation containing an uninstantiated

                          %  variable.

                          %

                          %  On exit, that variable is instantiated

                          %  to an approximate numeric solution.

                          %

                          %  The syntax of Left and Right is the same

                          %  as for expressions for the 'is' predicate.


                          solve(Left=Right) :-

                             free_in(Left=Right,X),

                             !,    /* accept only one solution of free_in */

                             define_dif(X,Left=Right),

                             solve_for(X).



                          % free_in(Term,Variable)

                          %

                          %  Variable occurs in Term and is uninstantiated.


                          free_in(X,X) :-                     % An atomic term

                             var(X).


                          free_in(Term,X) :-                  % A complex term

                             Term \== [[]],

                             Term =.. [_,Arg|Args],

                            (free_in(Arg,X) ; free_in(Args,X)).



                          % define_dif(X,Left=Right)

                          %

                          %  Defines a predicate to compute Left-Right

                          %  for the specified equation, given X.


                          define_dif(X,Left=Right) :-

                             abolish(dif,2),

                             assert((dif(X,Dif) :- Dif is Left-Right)).









                                                                                   9


                  Listing 3. Procedures to implement the secant method.



                  % solve_for(Variable)

                          %

                          %  Sets up arguments and calls solve_aux (below).


                          solve_for(Variable) :-

                             dif(1,Dif1),

                             solve_aux(Variable,1,Dif1,2,1).



                          % solve_aux(Variable,Guess1,Dif1,Guess2,Iteration)

                          %

                          %  Uses the secant method to find a value of

                          %  Variable that will make the 'dif' procedure

                          %  return a value very close to zero.

                          %

                          %  Arguments are:

                          %   Variable  -- Will contain result.

                          %   Guess1    -- Previous estimated value.

                          %   Dif1      -- What 'dif' gave with Guess1.

                          %   Guess2    -- A better estimate.

                          %   Iteration -- Count of tries taken.



                          solve_aux(cannot_solve,_,_,_,100) :-

                             !,

                             write('[Gave up at 100th iteration]'),nl,

                            fail.


                          solve_aux(Guess2,Guess1,_,Guess2,_) :-

                             close_enough(Guess1,Guess2),

                             !,

                             write('[Found a satisfactory solution]'),nl.


                          solve_aux(Variable,Guess1,Dif1,Guess2,Iteration) :-

                             write([Guess2]),nl,

                             dif(Guess2,Dif2),

                             Slope is (Dif2-Dif1) / (Guess2-Guess1),

                             Guess3 is Guess2 - (Dif2/Slope),

                             NewIteration is Iteration + 1,

                             solve_aux(Variable,Guess2,Dif2,Guess3,NewIteration).



                          % close_enough(X,Y)

                          %

                          %  True if X and Y are the same number to

                          %  within a factor of 0.0001.

                          %


                          close_enough(X,Y) :-









                                                                                                                    10


                     Quot is X / Y,

                             Quot > 0.9999,

                             Quot < 1.0001.

