The Kitchin Research Group

I have been exploring ways to use emacs-lisp to express scientific ideas. In this post, we explore a partial symbolic numeric solver in Emacs-lisp. This involves some syntactic developments to more clearly identify something we want to solve for and to then generate the code required to solve it.

In section The Newton solver you can find a simple implementation of a Newton solver in emacs-lisp. This function allows you to numerically solve equations that can be written in the form \(f(x) = 0\) for \(x\) given an initial guess. You write a function for \(f(x)\) and pass the function to the solver. This is a standard approach used in Python with fsolve, for example. Here is an example of solving a trivial problem: \(x - 4 = 0\) just to check that it works. We use a lambda function for \(f(x) = x - 4 = 0\). The answer is \(x=4\).

(newton-f (lambda (x) (- x 4)) 2)

That syntax is not too bad, but we have the whole lambda expression in there, and some repetition of what we want to solve for as an argument and in the function. It would be interesting if we could just have an expression that gets solved, e.g. (newton-f (- x? 4) 2) where x? indicates the thing to solve for.

We can do that! We can take an expression, flatten it and find the variable names that end with ?. We should check that there is only one, but for now we don't. Here is an example that does that. I use a nested expression here just to illustrate that the code finds the x? variable correctly.

(require 'dash)

(let ((body '((* (- x? 4) 1))))
  (loop for item in (-flatten body)
        if (and (symbolp item) (s-ends-with? "?" (symbol-name item)))
        collect item))

So, given an expression we can identify the unknown that should be the argument to a lambda function. So, we create a macro that takes that expression and constructs a function to solve it, then calls newton-f on it. The macro is syntactically useful here because we do not have to quote the expression. Here is that macro.

(defmacro solve (expression guess)
  `(newton-f
    (lambda ,(loop for item in (-flatten expression)
                   if (and (symbolp item) (s-ends-with? "?" (symbol-name item)))
                   collect item)
      ,expression)
    ,guess))

I call this a partial symbolic solver because we do some introspection symbolically to identify what to solve for, and then construct the code required to solve it. Here is that trivial example (x? - 4 = 0). It just shows we can have some nesting and it still works. I am not so thrilled with the initial guess, but this is an iterative solver, so you either need an initial guess, or a solution range.

(solve (* (- x? 4) 1) 3)

Here is what that expands into:

(macroexpand '(solve (* (- x? 4) 1) 3))

It expands into what we would have written in the first place. The benefit to us is less typing, and a simpler syntax. Both of those reduce the opportunity to make errors!

A more realistic problem might be: Reactant A flows into a continuously stirred tank reactor at a rate of \(F_{A0} = 1\) mol/min with a volumetric flow of \(v_0 = 1\) L/min.. The reactor achieves 50% conversion (\(X\)) of A to products. The reaction rate law is known to be \(-r_A = k C_A\) with \(k = 0.1\) 1/min. Estimate the volume of the reactor. If you have taken my class in reaction engineering, you know the following facts:

• The exit molar flow is defined by \(F_A = F_{A0} (1 - X)\)
• The exit concentration is \(C_A = F_A / v_0\)
• The mole balance is defined by \(0 = F_{A0} - F_A + r_A V\)

That is all we need; we can solve for \(V\) from the last equation. This is simple enough you might do the algebra to get: \(V = \frac{F_{A0} - F_A}{-r_A}\) which can be simply evaluated. We use our solver here and compare it to the evaluation.

Here is the solver:

(let* ((Fa0 1)
       (X 0.5)
       (Fa (* Fa0 (- 1 X)))
       (k 0.1)
       (v0 1)
       (Ca (/ Fa v0))
       (r (* k Ca))
       (ra (* r -1)))
  (solve (+ Fa0 (* Fa -1) (* ra V?)) 2))

It is pretty hard to imagine doing something like this in Python! It would probably involve parsing a string.

Here is the evaluation from our algebra:

(let* ((Fa0 1)
       (X 0.5)
       (Fa (* Fa0 (- 1 X)))
       (k 0.1)
       (v0 1)
       (Ca (/ Fa v0))
       (r (* k Ca))
       (ra (* r -1)))
  (/ (- Fa0 Fa) (* -1 ra)))

Within the tolerance specified in newton-f, these are the same.

This is just the tip of the iceberg though. You may have noticed that none of the variables in the let* had any descriptions. Sure, you could put some comments after them, but those are not really part of the code.

Also, we had to define the variables in advance of the expression. That is a limitation of how computers work, they cannot evaluate undefined variables. It constrains how we can express the idea. What if we could instead specify the equation first, then the data? That way we are clear what we are trying to do at a higher level, and fill in the details later. Suppose we wanted a syntax like the block below instead. Here we emphasize the equation we are solving first, and then define the variables and quantities used in the equation, and finally the guess that we use to find the solution.

(solve1
 (eqn (+ Fa0 (* -1 Fa) (* ra V?)))
 (data ((k 0.1 "rate constant 1/min")
        (Ca0 1.0 "feed concentration")
        (v0 1 "volumetric flow L/min")
        (Fa0 (* v0 Ca0) "Inlet molar flow")
        (X 0.5 "Desired conversion")
        (Fa (* Fa0 (- 1 X)) "Exit molar flow")
        (Ca (/ Fa v0) "exit concentration")
        (ra (* -1 k Ca) "rate in the reactor")))
 (guess 8))

That is achievable with the solve1 macro below! It too has some limitations, mostly the order of the data block still has to be correct, e.g. you cannot use a variable before it is defined. It would take some serious macro-fu to sort these so that everything is defined in the right order! Still, it allows you to express an executable idea in the order we defined. The strings in this syntax for documenting the variables are ignored, but they could be used in the macro to print useful information or something else you could imagine. You could also make them mandatory to encourage documentation.

(defmacro solve1 (&rest body)
  (let ((expression (second (assoc 'eqn body)))
        (data (loop for d in (second (assoc 'data body))
                    collect (list (first d) (second d))))
        (guess (second (assoc 'guess body))))
    `(let* ,data
       (newton-f
        (lambda ,(loop for item in (-flatten expression)
                       if (and (symbolp item) (s-ends-with? "?" (symbol-name item)))
                       collect item)
          ,expression)
        ,guess))))

To summarize, lisp macros allow us to rewrite the syntax of code before it is evaluated. This gives us the opportunity to inspect it, and generate new code, e.g. functions with arguments based on the contents of expressions, to save us typing. It also allows us to define ideas in a near arbitrary order that make sense to us, and then rearrange them so they make sense to the computer. See, for example, this post for an example of changing how functions are defined.

This seems to be heading in the domain specific language direction. I think it would be very helpful in engineering problem solving to build up tools like this. They could start out simple for new students to use. They never need to see the macro parts of this, just to learn how to use them for problem solving. These beginner tools would be limited in what they could do to minimize how much lisp is required to be learned so students can focus on the problem solving. Eventually they might outgrow them, and in the process transition to having the full lisp language at their disposal for problem solving.

;; See https://en.wikipedia.org/wiki/Newton%27s_method for the method

(defun newton-f (func x0)
  "Solve the equation FUNC(x)=0 using Newton's method.
X0 is an initial guess."
  (let* ((tolerance 1e-6)
         (x x0)
         (dx 1e-6)
         fx fpx)
    (while (> (abs (funcall func x)) tolerance)
      (setq fx (funcall func x)
            fpx (/ (- (funcall func (+ x dx)) (funcall func (- x dx))) (* 2 dx))
            x (- x (/ fx fpx))))
    x))

Copyright (C) 2017 by John Kitchin. See the License for information about copying.

org-mode source

Org-mode version = 9.0.5