Yesterday I wrote about Householder’s higher-order generalizations of Newton’s root finding method. For *n* at least 2, define

and iterate *H*_{n} to find a root of *f*(*x*). When *n* = 2, this is Newton’s method. In yesterday’s post I used Mathematica to find expressions for *H*_{3} and *H*_{4}, then used Mathematica’s `FortranForm[]`

function to export Python code. (Mathematica doesn’t have a function to export Python code per se, but the Fortran syntax was identical in this case.)

Aaron Muerer pointed out that it would have been easier to generate the Python code in Python using SymPy to do the calculus and `labdify()`

to generate the code. I hadn’t heard of `lambdify`

before, so I tried out his suggestion. The resulting code is nice and compact.

from sympy import diff, symbols, lambdify def f(x, a, b): return x**5 + a*x + b def H(x, a, b, n): x_, a_, b_ = x, a, b x, a, b = symbols('x a b') expr = diff(1/f(x,a,b), x, n-2) / \ diff(1/f(x,a,b), x, n-1) g = lambdify([x,a,b], expr) return x_ + (n-1)*g(x_, a_, b_)

This implements all the *H*_{n} at once. The previous post implemented three of the *H*_{n} separately.

The first couple lines of `H`

require a little explanation. I wanted to use the same names for the *numbers* that the function *H* takes and the *symbols* that SymPy operated on, so I saved the numbers to local variables.

This code is fine for a demo, but in production you’d want to generate the function `g`

once (for each *n*) and save the result rather than generating it on every call to `H`

.

Comments are closed.