Yesterday I wrote about Householder’s higher-order generalizations of Newton’s root finding method.
For n at least 2, define and iterate Hn 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 H3 and H4, 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 Hn at once.
The previous post implemented three of the Hn 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.
.