The Ultimate Algorithm for Solving Equations Numerically

数值方法解方程之终极算法 苏剑林 · 2010-04-04 · Mathematics equation series function numerical algorithm iteration taylor series

Ha, this title is a bit of a clickbait move on my part—maybe a touch exaggerated. I call it the "ultimate algorithm" mainly because it can be pushed to arbitrarily high precision and can, in principle at least, handle almost any nonlinear equation. Raising the precision just means adding more terms to an already-known iteration formula, rather than changing its form entirely—though of course as the precision goes up, so does the computational cost. The theoretical foundation here is, once again, the Taylor series.

Consider the equation $x=f(y)$. Given y, solving for x is easy, but given x, solving for y is not. Let's try expanding y as a Taylor series in x around $(x_0,y_0)$. The key is to work out the n-th derivative of y, $\frac{d^n y}{dx^n}$. Let us write $f^{(n)}(y)=\frac{d^n x}{dy^n}$, and we have

$$\frac{dy}{dx}=\frac{1}{(\frac{dx}{dy})}=f'(y)^{-1}$$more

and by applying the formula $y^{(n)}=y'\frac{dy^{(n-1)}}{dy}$, we can go on to obtain the second and third derivatives:

$$\begin{aligned}\frac{d^2 y}{dx^2}=\frac{d(f'(y)^{-1})}{dy}\cdot \frac{dy}{dx}=\frac{-f''(y)}{f'(y)^3} \\ \frac{d^3 y}{dx^3}=\frac{d(\frac{-f''(y)}{f'(y)^3})}{dy}\cdot \frac{dy}{dx}=\frac{3f''(y)^2}{f'(y)^5}-\frac{f'''(y)}{f'(y)^4}\end{aligned}$$

......

so that we can write

$$y=y_0+\frac{x-x_0}{f'(y_0)}-\frac{f''(y_0)}{f'(y_0)^3}\frac{(x-x_0)^2}{2!}+(\frac{3f''(y_0)^2}{f'(y_0)^5}-\frac{f'''(y_0)}{f'(y_0)^4})\frac{(x-x_0)^3}{3!}+...$$

and noting that $x_0=f(y_0)$, we have in effect derived the recursive formula for the root of the equation:

$$\begin{aligned}y_{n+1}=y_n+\frac{x-f(y_n)}{f'(y_n)}-\frac{f''(y_n)}{f'(y_n)^3}\frac{(x-f(y_n))^2}{2!}+ \\ (\frac{3f''(y_n)^2}{f'(y_n)^5}-\frac{f'''(y_n)}{f'(y_n)^4})\frac{(x-f(y_n))^3}{3!}+...\end{aligned}$$

In practice you don't need too many terms—just the first few will do (adding more terms causes the computation to blow up). If you keep only the first two terms, this reduces to the widely-used "Newton's method." Taking a few more terms speeds up convergence and improves stability.

Well, with this, Scientific Spaces has now introduced two "new" numerical algorithms for solving equations (the earlier one being the "tangent-parabola method"), and going forward I probably won't dig much further into numerical methods for algebraic equations. As BoJone sees it, this "ultimate algorithm" is more than enough for solving equations—paraphrasing Archimedes: "Give me an equation, and I shall pry out its root!"......