Solving Non-Linear Equations Using Computers

Non-linear equations are those that do not have all linear terms, i.e., each term may not be linear. There are many methods for solving these equations, including bisection, fixed point iteration, Newton's method, and the Secant method.

Solving an equation $f(x)$ is the equivalent of finding the roots of the equation, i.e., where $f(x) = 0$. For $f(x) = x^2 + 5x + 6$, for example, the roots are easy to solve as $f(x)$ factorises to $(x+2)(x+3)$, giving roots $-2,-3$.

For other equations, we may have complex roots. Take $f(x) = x^2 + 4x + 10$, for example. This has roots at $-2 \pm i\sqrt{6}$. Some equations do not permit solutions analytically, and thus must be computed using one of the various numerical techniques to find a sufficient approximation of the root.

Types of Equation

There are many different ways we can categorise equations. Any linear equation is known as a first order polynomial. Equations that pass through the $x$-axis multiple times are known as higher order polynomials. $\tan x - \sin x = 0$ would be an example of a transcendental equation.

Explicit vs Implicit

An explicit equation is where the equation is defined explicitly in terms of other variables and parameters, e.g., $y-mx+b = 0$ can be rearranged for $m$.

Implicit equations, on the other hand, do not permit us to isolate the variable on the left hand side of the equation.

Bisection

This is the process of working on a smaller and smaller problem, to get a root until such a point that we are very close to the root (given some $\epsilon = 0.001$, for example).

This is an iterative method. We choose a lower and upper guess for the root ($x_L$, $x_U$, respectively), such that the function changes signs over the interval.

This is easy to verify as $f(x_L)f(x_U) < 0$ if it does. We can then estimate the root at the midpoint of the guess $x_r = \frac{x_U + x_L}{2}$.

If the root is in the lower sub-interval, i.e., $f(x_L)f(x_r) < 0$, then we set the upper interval to $x_R$. It is trivial to see that if the root is in the upper sub-interval, then we go the other way.

If the $f(x_L)*f(x_r) = 0$, then we have found the root. This is because we converge on the root as $\lim_{x \to \infty}$, but may never equal 0 exactly. This is where the $\epsilon$ comes in handy.

Convergence Time

To find when we are within $\epsilon$ of the true solution $r$, then the error at each iteration $e_n$ is defined as:

\[ e_n = | r_n - r| \le \frac{x_R^n - x_L^n}{2} \le \frac{x_R^0 - x_L^0}{2^{n+1}} \]

This essentially says that the error is always less at each iteration. The latter part of the equation (after the last $\le$) shows that that term gives an upper bound for the time to converge.

For this error $e_n$ to be less than $n$, we assign another inequality:

\[ \begin{align} \frac{x_R^0 - x_L^0}{2^{n+1}} &\le \epsilon \\ x_R^0 - x_L^0 &\le 2^{n+1} \cdot \epsilon \\ \ln (x_R^0 - x_L^0) &\le \ln ( \epsilon \cdot 2^{n+1}) \\ \text{using rule } \ln(ab) &= \ln(a) + \ln(b):\\ &\le \ln(\epsilon) + \ln(2^{n+1})\\ &\le \ln(\epsilon) + (n+1)\ln(2)\\ &\le \ln(\epsilon) + n\ln(2) + \ln(2)\\ \text{using rule } \ln(ab) &= \ln(a) + \ln(b):\\ &\le (\ln(\epsilon) + \ln(2)) + n\ln(2) \\ &\le \ln(2\epsilon)+ n\ln(2) \\ \ln (x_R^0 - x_L^0) - \ln(2\epsilon) &\le n\ln(2) \\ \frac{\ln (x_R^0 - x_L^0) - \ln(2\epsilon)}{\ln(2)} &\le n \end{align} \]

This converges very slowly. This equation allows us to see the upper bound to the number of iterations that it will take for the solution to converge on an error tolerance, and this is useful before we do the computation to see if it will be time consuming to do.

Fixed Point Iteration

This is similar to bisection, but instead of continuing until we converge, we stop at a specific stop criterion. For this method, we rewrite $f(x) = 0$ in the form $x = g(x)$ (as $f(x) = 0 \Rightarrow f(x)+x=x)$.

Following this, we take a guess of the solution $x'$. If $g(x)$ is nice, then hopefully $g(x')$ is closer to $r$ than $x'$.

We start with $x_0$, then iterate $x_{k+1} = g(x_k)$ until we reach the stop criterion, which could be an error tolerance $\epsilon$, or a fixed number of iterations.

Example

Approximate $f(x) = x - \cos x$ to 4 digits of accuracy. We choose $g(x) = \cos x$, with $x = \cos x$. Start with $x_0 = 1$, then iterate $x_{k+1} = \cos x_k$.

After 25 iterations, we yield the same result to 4 digits of accuracy. Thus, our approximation is $x = 0.7391$.

Divergence

Sometimes, we will diverge when doing this, however, as the guess for $x$ may start to head off not in the direction of the root. The idea that the fixed point iteration method works is based on convergence, which in turn is based on $x_0$ and $g(x)$.

If we have continuous functions, and $r$ is the root $g(r) = r$, then the iteration $k$ is defined as $x_{k+1} = g(x_k)$. The error for this is therefore:

\[ \begin{align} e_{k+1} &= x_{k+1} - r\\ &= g(x_k) - g(r) \\ &= g'(\zeta)(x_k-r) \text{ such that } \zeta \in (x_k,r)\\ &= g'(\zeta) e_k\\ |e_{k+1}| &\le |g'(\zeta)||e_k| \end{align} \]

The error increases if $|g'| \ge 1$, leading to divergence, or decreases if $|g'| < 1$. From this, we can show that there is an interval of size $a$ around $r$ (this would be $[r-a, r+a] \land |g'(x)| <1$ for almost all $x$ in the interval).

For $g(x) = \cos x$, and thus $g'(x) = \sin x$, choosing $r=0.7391$ converges as $|g'(r)| < 1$. Here, $g'$ is simply the derivative of the function.

For a secondary example, $g(x) = 2^{-2x}(x-1)+x$, we can derive $g'(x) = -2e^{-2x}(x-1)+2e^{-2x}+1$. We can then take $|g'(1)| = 1 + e^{-2} > 1$. This shows divergence.

Convergence Iterations

If the iteration converges, then we can compute how many iterations are needed to get within an error tolerance $\epsilon$.

If we assume that $g'(x) \le m \in [r-a, r+a]$, then $|e_{k+1}| \le m|e_k|$, where $m$ is some constant. This is based on the above idea that if the gradient of the derivative decreases, then the error at the next iteration will be less than the error at the last iteration.

Using the $g'(x) \le m$ property, it is easy to see that $|e_1| \le m|e_0|$, and this gives:

\[ \begin{align} | e_1 | &\le m | e_0 |\\ | e_2 | &\le m | e_1 | \le m^2|e_0|\\ | e_k | &\le m^k| e_0 |\\ \end{align} \]

Warning

This part is unfinished, as I don't fully understand how the slide derives the maximum number of iterations required for convergence.

We could always replace $m$ with $g'(\zeta)$, but as we are working to an upper bound, we can just assume some constant $m$. This doesn't give the tightest possible bound, only some strict upper limit. We can tighten this bound using the following idea:

\[ \begin{align} |e_0| &= |r - x_0 | \text{ (the first error is the distance of first $x$ from $r$)}\\ &= |(r - x_1) + (x_0 - x_1)|\\ &\le |e_1| + |x_ - x_0|\\ &\le m|e_0| + |x_1 - x_0| \end{align} \]

Other than the derivation, all we need to understand is that the error $|e_k|$ is:

\[ | e_k | \le \frac{m^k}{1-m} | x_1 - x_0 | \]

so to get in some distance $\epsilon$, we need $\epsilon \ge |e_k|$:

\[ \frac{m^k}{1-m} | x_1 - x_0 | \le \epsilon \]

Then we can obtain $k$ by rearranging:

\[ \begin{align} \frac{m^k}{1-m} |x_1 - x_0| &\le \epsilon\\ m^k|x_1 - x_0| &\le \epsilon(1-m) \\ \ln(m^k |x_1 - x_0|) &\le \ln(\epsilon(1-m)) \\ \ln(m^k) + \ln(|x_1 - x_0|) & \le \ln( \epsilon ( 1 -m))\\ \ln(m^k) &\le \ln( \epsilon ( 1 -m)) - \ln(|x_1 - x_0|)\\ k \ln(m) &\le \ln( \epsilon ( 1 -m)) - \ln(|x_1 - x_0|)\\ k &\le \frac{\ln( \epsilon ( 1 -m)) - \ln(|x_1 - x_0|)}{\ln(m)}\\ \end{align} \]

Taking the following example, we can see how many iterations it will take to converge for $x_0 = 1$:

\[ \begin{align} f(x) &= \cos (x) - x \\ g(x) &= \cos (x) \\ g'(x) &= -\sin(x) \\ x_0 &= 1\\ \\ \text{check }|g'(x)| &\le 1\\ &= |-\sin(1)| = 0.8415\dots = m\\ \\ x_1 &= \cos (x_0)\\ &= \cos (1)\\ &= 0.5403\dots\\ \\ \text{so for }\epsilon &\le 10^{-5}\\ k &= \frac{\ln(10^{-5})(1-0.8415\dots)) - \ln(|0.5403\dots - 1|)}{\ln m}\\ &\approx 73 \end{align} \]

Newton's Method

Newton's Method considers $f(x) = 0$, trying to get the root of a function. It starts with an initial guess $x_0$. To get $x_{k+1}$ from $x_k$, we treat $f$ as a linear function, and determine the root of that function. It is essentially an optimal version of FPI.

The mathematical way of doing this is to linearise $f(x)$ at $x_k$, giving $y(x)$:

\[ y(x) = \underbrace{f(x_k)}_{y\text{-intercept}} + \underbrace{(x-x_k)}_{\text{shift left}}\underbrace{f'(x_k)}_{\text{gradient at }x_k} \]

Essentially, convert this into a $y = mx + c$, where the $c=f(x_k)$, and the $mx$ is simply the slope of the line.

The root of this function is then when $y(x) = 0$, which is at $x = x_k -f(x_k)/f'(x_k)$. We use the root of $y$ as $x_{k+1}$ and then iterate, creating the following iterative formula:

\[ x_{k+1} = x_k - f(x_k)/f'(x_k) \]

We can also show that Newton's method is a fixed point iteration, as $g(x) = x - f/f'$. To help with this, we introduce the notion of the $b(x)$ function as a 'relaxation factor'. This is chosen such that we can rewrite $f(x) = 0$ as

\[ \begin{align} -b(x)f(x) &= 0\\ \Rightarrow g(x) &= -b(x)f(x) = x \end{align} \]

We can then take the derivative of $g(x)$ with respect to $x$, yielding:

\[ g'(x) = 1 - b'(x)f(x) - f'(x)b(x) \]

From earlier notes, we know that the lower $g'(x)$ is, the faster the function will converge, and if $g'(x)> 1$, then the function won't converge at all. The function will converge fastest if we select $b(x) = 1/f'(x)$.

Convergence Iterations

This is largely the same as the analysis for FPI, with $r$ being the true root and $e_{k+1} = x_{k+1} - r$, which is $g(x_k) - g(r)$. We can take the Taylor expansion of this function, stopping after the second derivative. Bear in mind the Taylor series:

\[ f(a) + \frac {f'(a)}{1!}(x-a) + \frac{f''(a)}{2!}(x-a)^2 \]

Substituting $g$ for $f$ and $x_k$ for $a$, and $r$ for $x$: $$ g(x_k) = g(r) + \frac {g'(r)}{1}(x_k-r) + \frac{g''(r)}{2}(x_k-r)^2 $$

The gradient $g'(r) = 0$ allows us to simplify to $g(x_k) = g(r) + \frac{g''(\zeta)}{2} (x_k - r)^2 | \zeta \in (x_k, r)$, giving:

\[ e_{k+1} = \frac{1}{2}e_k^2|g''(\zeta)| \]

Using the Newton method for an example $f(x) = x^2 - a$, we see that

\[ \begin{align} x_{k+1} &= x_k + \frac{f(x_k)}{f'(x_k)} \\ &= \frac{x_k}{2} + \frac{a}{2x_k} \end{align} \]

So, for an example $a=3, x_0 = 1.7$, we have $x_0 = 1.7$, $x_1 = 1.7324$, $x_2 = 1.7321$, $x_3 = 1.7321$.

Problems

The method converges slowly for some functions (e.g., if there is an inflection point near a root). Local minima and maxima can also lead to oscillation.

If a function has multiple roots the method also struggles to deal with it. If this is the case, then we can use a new function $u = f'/f$, which has roots at the same $x$ as $f(x)$. In this instance, we can use Newton with the modified function.

Secant Method

If the function is complicated, then it may be hard to find $f'(x)$. We can use the Newton scheme and approximate $f'(x)$ using the secant:

\[ \begin{align} f'(x_k) &\approx \frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}}\\ \text{substitute into Newton:}\\ x_{k+1} &\approx x_k - \frac{x_k - x_{k-1}}{f(x_k) - f(x_{k-1})}f(x_k)\\ \end{align} \]

This provides us with a way to approximate the Newton method for a complicated function. We can also show that the error $e_{k+1} \le Ce_k^\alpha$, where $\alpha = \frac{1}{2}(1+\sqrt{5})$, but this was not derived in the lecture. It is a bit slower than Newton, but it means that we can still use something similar to Newton for approximating the roots.

Practical Approaches to Root Finding

We have discussed the benefits and problems of each of the methods above. Typically, a program will use 5-6 iterations of a bisection to find a good initial guess of $x_0$. Once this has been found, then we can use Newton's method for the next 3-4 iterations to refine that guess, as Newton is quite good when $x_0$ close to the real solution.

Systems of Equations

We have so far discussed solving one equation with a single $x$. We typically have this in a higher dimension. The system would be $\vec{F}(\vec{x}) = 0$, where $\vec{F} = (f_1, \dots, f_n), \vec{x} = (x_1, \dots, x_n)$

We can rewrite these in terms of $\vec{G}$ for fixed-point iterations and choose $\vec{x_0}$. Using this, we then see that Newton's method generalises to:

\[ \vec{x_{k+1}} = \vec{x_k} - (DF(\vec{x_k}))^{-1} \vec{F}(x_k) \]

But hang on! What's this $DF(\vec{x_k})$? It's the matrix inverse of the total derivative or Jacobian matrix. Essentially each row of the matrix is each $f_x$ and each column works through the variables $x_1, \dots, x_n$

We can now solve equations with multiple unknowns! Consider:

\[ \begin{align} x^2y + y &= 0\\ y \sin(x) + xy^2 &=0\\ \\ \vec{x} &= (x,y)\\ F(\vec{x}) &= (x^2y+y, y\sin(x) + xy^2)\\ \\ DF(\vec{x}) &= \begin{pmatrix} 2xy & x^2 + 1\\ y\cos(x) + y^2 & \sin(x) + 2xy\\ \end{pmatrix} \end{align} \]

I won't attempt to write down the inverse of $DF$, just know that it is possible to do and that it is horrible even for a matrix rank 2.