Nonlinear equation I

  • Find roots of $f(x)$. That is, find the value of $x$ so that $f(x)=0$.
  • All nonlinear algebraic problems can be written in the form $f(x)=0$. Just move all terms to one side of the equation.
  • Numerical methods for nonlinear problems proceed by setting an initial guess for the solution, and iterating to improve the guess until some desired error tolerance is achieved.
  • Open and Closed domain methods.
  • Here:
    • discuss closed domain methods
    • one equation in one unknown

Bound the solution

  • Bracket the root:
    • Choose two guesses for the solution: $x_1$ and $x_2$. If $f(x_1)\cdot f(x_2)<0$ then the two guesses bracket the root.
    • This implies that the root is between the two guesses.

Procedures

  • Graph it.
    • Visual picture can be very helpful.
    • Provides an intial guess of the root.
    • Indicates the function behavior, and possible problem areas.
    • May not be practical, however.
      • Cost of function evaluation may be excessive.
  • Do a simple incremental search.
    • Simply set a small $\Delta x$ and search the domain for $f(x)=0$.
    • Obviously not very practical.
  • Past experience
    • Output of one problem may be a good guess for the input of another problem.
      • This especially true in cases where we do multiple solves, which is quite common.
    • Use intuition:
      • species mass fractions should be between 0 and 1.
      • most temperatures in K should be in reasonable ranges.
  • Solve a simpler problem to get a guess for the harder problem.
    • For example, if solving for a nonideal gas, get an initial guess using an ideal gas.
    • Another example, if solving for temperature with a variable heat capacity, find a temperature guess using a constant heat capacity, which results in an easy linear solve for the guess.
    • In combustion, we might solve using a simple 1-step reaction mechanism, instead of a more complex 400 step mechanism.

Once you have a bracket refine the root

  • Reduce the size of the interval, while maintaining the bracket.
  • Two methods: bisection, and regula falsi (false position)
    • These are robust (they work!), but they are not overly fast.

Bisection

  • Guess two points that bracket the root: $a$, $b$
  • Check for $f(a)=0$ or $f(b)=0$
  • Bracket tests:
    • $f(a)f(b)<0$
    • ($f(a)>0$ and $f(b)<0$) or ($f(a)<0$ and $f(b)>0$)
  • Refine: $c = (a+b)/2$
  • Select new bracket:
    • if $f(a)f(c) < 0 \rightarrow b=c$
    • else $a=c$

Error

  • The error is always bounded by $\epsilon_n\le|b-a|$
  • $\epsilon_{n+1} = \epsilon_n/2$ $\rightarrow$ linear convergence.
  • $\epsilon_0$, $\epsilon_1=\epsilon_0/2$, $\epsilon_2 = \epsilon_1/2 = \epsilon_0/4 = \epsilon_0/2^2$
  • $\rightarrow \epsilon_n = \epsilon_0/2^n$.
  • Hence: $$n = \log_2\left(\frac{\epsilon_0}{\epsilon_n}\right)$$
  • That is, to reduce the error from $\epsilon_0 \le |a-b|$ to some desired $\epsilon_n$ requires $n=\log_2(\epsilon_0/\epsilon_n)$ iterations.

Regula Falsi

  • Rather than bisect the interval to find $c$, draw an intersecting line between $f(a)$ and $f(b)$ as an approximation to the function.

  • Solve the root of this linear approximation:

    • that is, where the line intersects the x-axis where $y=f(x)$ is zero.
    • Equation for the line: $$f_l(x) = f(a) + \frac{f(b)-f(a)}{b-a}(x-a).$$
    • Then let $f_l(c)=0$ and solve for $c$ to get: $$ c= a - f(a)\frac{b-a}{f(b)-f(a)}.$$
  • To retain the bracket, replace $a$ or $b$ with $c$ as for bisection.

  • Robust
  • Usually faster than bisection
  • No error bound though.
  • May get superlinear convergence
    • Keeping the old versus new function evaluation…
  • Consider the following case though (slow):

Convergence

  • When to stop?
    • $|b-a| < \epsilon$
      • absolute error
      • works for $x=1$
      • not so good for $x=10^{20}$
    • $|b-a|/|c| < \epsilon$
      • relative error
    • $|f(c)|< \epsilon$
      • error in function versus error in root.

Consider:

  • On the left, the bracket is narrow, but the function is not zero.
  • On the right, the function is near zero, but the bracket is not narrow.
  • Do both $|b-a|/|c| < \epsilon_1$ and $|f(c)|<\epsilon_2$.

Bisection example

Find root of $f(x)=(x-2)^2 - 1$.

import numpy as np
def bisection(f, a, b, tol=1E-5, maxit=10000):
    
    fa = f(a)
    fb = f(b)

    if fa==0: return a, 0
    if fb==0: return b, 0

    if fa*fb > 0.0:
        print("a, b don't bracket the root, try again")
        return np.nan, -1

    for k in range(1,maxit+1):

        c = 0.5*(a+b)
        fc = f(c)
        if fc==0: return c, k

        if fa*fc <= 0: 
            b = c
        else:
            a = c

        if np.abs(fc) <= tol or np.abs((b-a)/c) <= tol:
            return c, k

    print(f"no convergence in {maxit} iterations")
    return c, maxit
        
#--------------------        

def func1(x):
    return (x-2)**2 - 1

#--------------------        

x, nit = bisection(func1, 2., 3.5, tol=1E-8, maxit=100)

print(f"Solution x = {x} found in {nit} iterations")
print(f"f(x) = {func1(x)}")
Solution x = 3.0000000074505806 found in 26 iterations
f(x) = 1.4901161193847656e-08

Compare with fsolve

from scipy.optimize import fsolve

xguess = 2
x = fsolve(func1, xguess)[0]

print(f"Solution x = {x}")
Solution x = 3.0