# The Solution of Nonlinear Equations

<a rel="license" href="http://creativecommons.org/licenses/by/4.0/"></a><br />This notebook by Xiaozhou Li is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by/4.0/">Creative Commons Attribution 4.0 International License</a>.  
All code examples are also licensed under the [MIT license](http://opensource.org/licenses/MIT).

One of the most frequently occurring problems in scientific work is to find the roots of equations of the form
\begin{equation}
    f(x) = 0,
\end{equation}
i.e., zeros of the function $f(x)$.

## The Existence of Root
The first step to solving an equation is to verify that a root exists.

### Theorem
   Let $f$ be a continuous function on $[a, b]$, satisfying $f(a)f(b) < 0$. Then $f$ has a root between $a$ and $b$, that is, there exists a number $x^{*}$ satisfying $a < x^{*} < b$ and $f(x^{*}) = 0$.
   
__Notes:__ 
* This theorem is just the simplest one for the existence of roots. 

## The Bisection Method
Let us see one of the simplest root finding methods: the bisecion method.

### Bisection Method 
**Algorithm:**
```python
    Given initial interval [a, b] such that f (a)f (b) < 0 
    while (b − a)/2 > TOL:
        c = (a + b)/2
    if f (c) = 0:
        stop 
    if f(a)f(c) < 0:
        b=c 
    else:
        a=c
    The final interval [a,b] contains a root. 
    The approximate root is (a + b)/2.
```

**Python implementation:**


```python
# environment setting, before any codes
import numpy as np
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import matplotlib.pyplot as plt
plt.rcParams['axes.unicode_minus']=False
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
from IPython.display import clear_output, display
```

### The Python Code
Here, we give the code of bisection method as follows:
```python
    def bisect(f, a, b, tol):
```
* $f$ is the function $f(x)$ 
* $a$ and $b$ are for the interval $[a, b]$
* $tol$ is the tolerence error for the numerical solution. 


```python
def bisect(f, a, b, tol):
    if (np.sign(f(a)*f(b)) >= 0):      # wrong input
        print ('$f(a)f(b) <0$ not satisfied')
        return
    fa = f(a) 
    fb = f(b)
    while ((b - a)/2) > tol:
        c = (a + b)/2;
        fc = f(c)
        if fc == 0:                 # c is a solution, done
            return c 
        if (np.sign(fc*fa) < 0):       # new interval [a, c]
            b = c; 
            fb = fc
        else:                       # new interval [c, b]
            a = c
            fa = fc
    return (a + b)/2                #new midpoint is best estimate
```

### Demostration
**Demonstration version:**


```python
def bisect_iters(f, a, b, iters):
    fa = f(a) 
    fb = f(b)
    
    x = np.zeros(iters+1)
    for i in range(iters):
        c = (a + b)/2
        x[i] = c
        fc = f(c)
        if fc == 0:                 # c is a solution, done
            return c 
        if (np.sign(fc*fa) < 0):       # new interval [a, c]
            b = c; 
            fb = fc
        else:                       # new interval [c, b]
            a = c
            fa = fc
    x[iters] = (a + b)/2       
    return x              #new midpoint is best estimate
```

**Exmaple**
    Use Bisection Method to find the root of $f(x) = e^{-x} - \sin\left(\frac{\pi x}{2}\right)$ on $[0,1]$.


```python
def fun(x):
    return np.exp(-x) - np.sin(np.pi/2*x)
```


```python
def Iteration_demo(f, a, b, iters):
    x = np.linspace(a, b, 1000)
    plt.plot(x, f(x), 'k', linewidth=3)
    plt.plot(x, 0*x, 'k-.', linewidth=3)
    x = bisect_iters(f, a, b, iters)
    x_ref = bisect(f, a, b, 1.e-14)
    
    #print (x)
    plt.xlim(a,b)
    plt.ylim(min(f(a),f(b)),max(f(a),f(b)))
    #plt.ylim(1,2)
    for i in range(iters+1):
        plt.plot(x, f(x), 'r*', markersize=12, linewidth=3)
    
    print ('Reference Solution: ', x_ref)
    print ('no.    solution    error bound    error')
    for i in range(iters+1):
        print ("%3d  %12.10f   %7.2e    %7.2e" % \
                   (i, x[i], (b-a)/2**(i+1), np.abs(x_ref - x[i])))
    
        
a = 0
b = 1
Iteration_demo(fun, a, b, 10)
```


```python
w = interactive(Iteration_demo, f=fixed(fun), a=fixed(0), b=fixed(1), iters=widgets.IntSlider(min=0,max=50,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=50), Output()), _dom_classes=('widget-intera…


**Exmaple**
    Use Bisection Method to find the root of $f(x) = x^3 + x - 1$ on $[0,1]$.


```python
def fun2(x):
    return x**3 + x - 1

w = interactive(Iteration_demo, f=fixed(fun2), a=fixed(0), b=fixed(1), iters=widgets.IntSlider(min=0,max=20,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=20), Output()), _dom_classes=('widget-intera…


**Exmaple**
    Use Bisection Method to find the root of __linear equation__ $f(x) = x - 1/3$ on $[0,1]$.  What do you think of the performance?


```python
def fun3(x):
    return x - 0.50000000001

w = interactive(Iteration_demo, f=fixed(fun3), a=fixed(0), b=fixed(1), iters=widgets.IntSlider(min=0,max=20,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=20), Output()), _dom_classes=('widget-intera…


## Fixed-Point Iteration (FPI)

### FPI Algorithm
First, one can rewrite the equation $f(x) = 0$ as a fixed point problem $g(x) = x$, Fixed-Point Iteration proceeds by starting with an initial guess $x_0$ and iterating the function $g$.
\begin{align}
    x_0 & = \text{initial guess} \\
    x_{n+1} & = g(x_n) \,\,\text{for}\,\,n = 0, 1,\ldots,
\end{align}

### The Python Code
The code of FPI is very easy, in each iteration, just do a function evaluation.


```python
def fixed_iters(g, x0, iters):
    x = np.zeros(iters+1)
    x[0] = x0
    for i in range(iters):
        x[i+1] = g(x[i])
    return x

def g(x):
    return np.cos(x)

fixed_iters(g, 1., 20)
```




    array([1.        , 0.54030231, 0.85755322, 0.65428979, 0.79348036,
           0.70136877, 0.76395968, 0.72210243, 0.75041776, 0.73140404,
           0.74423735, 0.73560474, 0.74142509, 0.73750689, 0.74014734,
           0.7383692 , 0.7395672 , 0.73876032, 0.73930389, 0.73893776,
           0.7391844 ])



### Demostration
**Demonstration version:**


```python
def Iteration_demo2(f, a, b, g, x0, iters):
    x = np.linspace(a, b, 1000)
    plt.plot(x, f(x), 'k', linewidth=3)
    plt.plot(x, 0*x, 'k-.', linewidth=3)
    x = fixed_iters(g, x0, iters)
    x_ref = bisect(f, a+1.e-8, b, 1.e-14)
    
    #print (x)
    plt.xlim(a,b)
    plt.ylim(min(f(a),f(b)),max(f(a),f(b)))
    for i in range(iters+1):
        plt.plot(x, f(x), 'r*', markersize=12, linewidth=3)
    
    print ('Reference Solution: ', x_ref)
    print ('no.    solution    |x_{n+1} - x_{n}|    error')
    for i in range(iters+1):
        if (i < iters):
            print ("%3d  %14.12f      %7.2e       %7.2e" % \
                   (i, x[i], np.abs(x[i+1]-x[i]), np.abs(x_ref - x[i])))
        else:
            print ("%3d  %14.12f      %7.2e       %7.2e" % \
                   (i, x[i], np.abs(g(x[i])-x[i]), np.abs(x_ref - x[i])))
```

**Example** Use Fixed Point Iteration to find a root of $f(x) = x^3 + 4x^2 - 10$ on $[1, 2]$ with $g_1(x) = \sqrt{\frac{10}{x} - 4x}$ and $g_2(x) = \sqrt{\frac{10}{x+4}}$, the initial guess $x_0 = 1.5$.


```python
def f(x):
    return x**3 + 4*x**2 - 10

a = 1; b = 2;

def g1(x):
    return np.sqrt(10/x - 4*x)

def g2(x):
    return np.sqrt(10/(x+4))

x0 = 1.5

w = interactive(Iteration_demo2, f=fixed(f), a=fixed(a), b=fixed(b), g=fixed(g2), x0 = fixed(x0), iters=widgets.IntSlider(min=0,max=20,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=20), Output()), _dom_classes=('widget-intera…


__Better than bisection?__
* Solving the above example by using bisection method.


```python
w = interactive(Iteration_demo, f=fixed(f), a=fixed(1), b=fixed(2), iters=widgets.IntSlider(min=0,max=50,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=50), Output()), _dom_classes=('widget-intera…


**Example** Find the fixed points of $g(x) = 2.8x - x^2$


```python
def f(x):
    return 1.8*x - x**2

a = 0; b = 2;

def g(x):
    return 2.8*x - x**2

x0 = -1

w = interactive(Iteration_demo2, f=fixed(f), a=fixed(a), b=fixed(b), g=fixed(g), x0 = fixed(x0), iters=widgets.IntSlider(min=0,max=40,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=40), Output()), _dom_classes=('widget-intera…


**Example** Use Fixed Point Iteration to find the positive root of $f(x) = x^2 - 3$ with $g_1(x) = x - \frac{1}{4}(x^2-3)$ and $g_2(x) = \frac{1}{2}\left(x + \frac{3}{x}\right)$.


```python
def f(x):
    return x**2 - 3

a = 0; b = 2

def g1(x):
    return x - 0.25*(x**2-3)

def g2(x):
    return 0.5*(x + 3/x)

x0 = 2

w = interactive(Iteration_demo2, f=fixed(f), a=fixed(a), b=fixed(b), g=fixed(g1), x0 = fixed(x0), iters=widgets.IntSlider(min=0,max=40,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=40), Output()), _dom_classes=('widget-intera…


## Newton's Method
### Newton's Method 
Newton's method is a particular kind o af fixed point iteration, with $g(x) = x - \frac{f(x)}{f'(x)}$.
\begin{align}
    x_0 & = \text{initial guess} \\
    x_{n+1} & = x_n - \frac{f(x_n)}{f'(x_{n})} \,\,\text{for}\,\,n = 0, 1,\ldots,
\end{align}

### The Python Code
Similar as the code of FPI, in each iteration, the Newton's method just need do a function evaluation and compute its derivative.


```python
def Newton_iters(f, f_der, x0, iters):
    x = np.zeros(iters+1)
    x[0] = x0
    for i in range(iters):
        x[i+1] = x[i] - f(x[i])/f_der(x[i])
    return x

def f(x):
    return x**2 - 3

def f_der(x):
    return 2*x

Newton_iters(f, f_der, 1., 20)
```




    array([1.        , 2.        , 1.75      , 1.73214286, 1.73205081,
           1.73205081, 1.73205081, 1.73205081, 1.73205081, 1.73205081,
           1.73205081, 1.73205081, 1.73205081, 1.73205081, 1.73205081,
           1.73205081, 1.73205081, 1.73205081, 1.73205081, 1.73205081,
           1.73205081])



**Example** Find the positive root of $f(x) = x^2 - a$ by Newton's Methods

The iterations is 
\begin{equation}
    x_{n+1} = \left(x_n + \frac{a}{x_n}\right)/2 
\end{equation}


```python
def Iteration_demo3(f, f_der, a, b, x0, x_star, iters):
    x = np.linspace(a, b, 1000)
    plt.plot(x, f(x), 'k', linewidth=3)
    plt.plot(x, 0*x, 'k-.', linewidth=3)
    x = Newton_iters(f, f_der, x0, iters)
    #x_ref = bisect(f, a+1.e-8, b, 1.e-14)
    
    #print (x)
    plt.xlim(a,b)
    #plt.ylim(min(f(a),f(b)),max(f(a),f(b)))
    for i in range(iters+1):
        plt.plot(x, f(x), 'r*', markersize=12, linewidth=3)
    
    print ('Reference Solution: ', x_star)
    print ('no.    solution    |x_{n} - x_{*}|    e_{n}/e_{n-1}^2   e_{n}/e_{n-1}')
    for i in range(iters+1):
        if (i < 1):
            print ("%3d  %14.12f      %7.2e" % \
                   (i, x[i], np.abs(x_star - x[i])))
        else:
            print ("%3d  %14.12f      %7.2e           %6.4f          %6.4f" % \
                   (i, x[i], np.abs(x_star - x[i]), np.abs(x_star - x[i])/(x_star - x[i-1])**2, \
                   np.abs(x_star - x[i])/np.abs(x_star - x[i-1])))
```


```python
def f(x): 
    return x**2 - 2

def f_der(x):
    return 2*x

x0 = -3
x_star = np.sqrt(2)
w = interactive(Iteration_demo3, f=fixed(f), f_der=fixed(f_der), a=fixed(a), b=fixed(b), x0 = fixed(x0), \
                x_star = fixed(x_star), iters=widgets.IntSlider(min=0,max=20,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=20), Output()), _dom_classes=('widget-intera…


**Example** Find the root of $f(x) = x^3 + 10x - 20$ around $x_0 = 1.5$ by Newton's Methods

The iterations is 
\begin{equation}
    x_{n+1} = x_n - \frac{x^3 + 10x - 20}{3x^2 + 10}
\end{equation}


```python
def f(x): 
    return x**3 + 10*x - 20

def f_der(x):
    return 3*x**2 + 10

x0 = 100
a = 1 
b = 2 
x_star = bisect(f, a, b, 1.e-14)
w = interactive(Iteration_demo3, f=fixed(f), f_der=fixed(f_der), a=fixed(a), b=fixed(b), x0 = fixed(x0), \
                x_star = fixed(x_star), iters=widgets.IntSlider(min=0,max=10,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=10), Output()), _dom_classes=('widget-intera…


**Example** Find the root of $f(x) = x^3 - 3x + 2$ around $x_0 = 1.1$ by Newton's Methods

The iterations is 
\begin{equation}
    x_{n+1} = x_n - \frac{x^3 - 3x + 2}{3x^2 - 3}
\end{equation}


```python
def f(x): 
    return x**3 - 3*x +2

def f_der(x):
    return (3*x**2 - 3)/1.

x0 = -1.5
a = -3. 
b = 3. 
x_star = -2.
w = interactive(Iteration_demo3, f=fixed(f), f_der=fixed(f_der), a=fixed(a), b=fixed(b), x0 = fixed(x0), \
                x_star = fixed(x_star), iters=widgets.IntSlider(min=0,max=30,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=30), Output()), _dom_classes=('widget-intera…


**Modified Newton Method**  
Find the root of $f(x) = x^3 - 3x + 2$ around $x_0 = 1.1$ by Modified Newton Methods

The iterations is 
\begin{equation}
    x_{n+1} = x_n - \frac{x^3 - 3x + 2}{x^2 - 1}
\end{equation}


```python
def f(x): 
    return x**3 - 3*x +2

def f_der(x):
    return (3*x**2 - 3)/2.

x0 = 1.1
a = -3. 
b = 3. 
x_star = 1.
w = interactive(Iteration_demo3, f=fixed(f), f_der=fixed(f_der), a=fixed(a), b=fixed(b), x0 = fixed(x0), \
                x_star = fixed(x_star), iters=widgets.IntSlider(min=0,max=30,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=30), Output()), _dom_classes=('widget-intera…


### Secant Method
One drawback of Newton's method is that it needs the function $f(x)$ and its derivative function $f'(x)$.  What if you do not have the derivative function $f'(x)$?  One way to avoid this issue is to use the so-called secant method.

__Secant Method:__
\begin{align}
    x_0, x_1 & = \text{initial guesses} \\
    x_{n+1} & = x_n - \frac{f(x_n)(x_n - x_{n-1}}{f(x_{n})-f(x_{n-1})} \,\,\text{for}\,\,n = 1, 2,\ldots,
\end{align}

__The Python Code__


```python
def secant_iters(f, x0, x1, iters):
    x = np.zeros(iters+1)
    x[0] = x0
    if (iters >= 1):
        x[1] = x1
        for i in range(1,iters):
            x[i+1] = x[i] - f(x[i])*(x[i]-x[i-1])/(f(x[i]) - f(x[i-1]))
    return x

def f(x):
    return x**2 - 3

secant_iters(f, 1., 2., 10)
```




    array([1.        , 2.        , 1.66666667, 1.72727273, 1.73214286,
           1.73205068, 1.73205081, 1.73205081, 1.73205081, 1.73205081,
           1.73205081])




```python
def Iteration_demo4(f, a, b, x0, x1, x_star, iters):
    x = np.linspace(a, b, 1000)
    plt.plot(x, f(x), 'k', linewidth=3)
    plt.plot(x, 0*x, 'k-.', linewidth=3)
    x = secant_iters(f, x0, x1, iters)
    #x_ref = bisect(f, a+1.e-8, b, 1.e-14)
    
    #print (x)
    plt.xlim(a,b)
    #plt.ylim(min(f(a),f(b)),max(f(a),f(b)))
    for i in range(iters+1):
        plt.plot(x, f(x), 'r*', markersize=12, linewidth=3)
    
    print ('Reference Solution: ', x_star)
    print ('no.    solution    |x_{n} - x^{*}|    e_{n}/e_{n-1}^1.618   e_{n}/e_{n-1}')
    for i in range(iters+1):
        if (i < 1):
            print ("%3d  %14.12f      %7.2e" % \
                   (i, x[i], np.abs(x_star - x[i])))
        else:
            print ("%3d  %14.12f      %7.2e           %6.4f          %6.4f" % \
                   (i, x[i], np.abs(x_star - x[i]), np.abs(x_star - x[i])/np.abs(x_star - x[i-1])**1.618, \
                   np.abs(x_star - x[i])/np.abs(x_star - x[i-1])))
```


```python
def f(x): 
    return x**3 - 3*x +2

x0 = -1.5
x1 = -2.5
a = -3. 
b = 3. 
x_star = -2

w = interactive(Iteration_demo4, f=fixed(f), a=fixed(a), b=fixed(b), x0 = fixed(x0), x1 = fixed(x1), \
                x_star = fixed(x_star), iters=widgets.IntSlider(min=0,max=20,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=20), Output()), _dom_classes=('widget-intera…


## Assignment
### Assignment 1:
For equation $x^2−3x+2−e^x =0$
1) Fixed Point Iteration: 
    \begin{equation*}
        x_{n+1} = g(x_n)
    \end{equation*}
     with $g(x) = (x^2 + 2^x + 2 - e^x)/5$

2) Newton's Methods: 
\begin{equation*}
x_{n+1} = x_n - \frac{x^2−3x+2−e^x}{2x -3 - e^x}
\end{equation*}


```python
def f(x):
    return x**2 - 3*x +2 - np.exp(x)

def g(x):
    return x + 0.2*(x**2 - 3*x +2 - np.exp(x))

a = 0; b = 1;
x0 = 0.5

w = interactive(Iteration_demo2, f=fixed(f), a=fixed(a), b=fixed(b), g=fixed(g), x0 = fixed(x0), iters=widgets.IntSlider(min=0,max=20,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=20), Output()), _dom_classes=('widget-intera…



```python
def f(x): 
    return x**2 - 3*x +2 - np.exp(x)

def f_der(x):
    return 2*x - 3 - np.exp(x)

x0 = 0.5
a = 0. 
b = 1. 
x_star = bisect(f, a, b, 1.e-14)
w = interactive(Iteration_demo3, f=fixed(f), f_der=fixed(f_der), a=fixed(a), b=fixed(b), x0 = fixed(x0), \
                x_star = fixed(x_star), iters=widgets.IntSlider(min=0,max=30,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=30), Output()), _dom_classes=('widget-intera…


### Assignment 2: 
Design a fixed point iteration method to solve equation $x^3 +2x^2 +10x −20 = 0$


```python
def f(x):
    return x**3 + 2*x**2 + 10*x - 20

def g(x):
    return x - 0.02*(x**3 + 2*x**2 + 10*x -20)

a = 1; b = 2;
x0 = 1.5

w = interactive(Iteration_demo2, f=fixed(f), a=fixed(a), b=fixed(b), g=fixed(g), x0 = fixed(x0), iters=widgets.IntSlider(min=0,max=30,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=30), Output()), _dom_classes=('widget-intera…



```python
def f(x): 
    return x**3 + 2*x**2 + 10*x - 20

def f_der(x):
    return 3*x**2 + 4*x + 10

x0 = 1.5
a = 1. 
b = 2. 
x_star = bisect(f, a, b, 1.e-14)
w = interactive(Iteration_demo3, f=fixed(f), f_der=fixed(f_der), a=fixed(a), b=fixed(b), x0 = fixed(x0), \
                x_star = fixed(x_star), iters=widgets.IntSlider(min=0,max=30,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=30), Output()), _dom_classes=('widget-intera…


### Assignment 3
Consider Fixed-Point Iteration applied to $g(x) = 1 − 5x + \frac{15}{2}x^2 − \frac{5}{2}x^3$ . (a) Show that 
$1 − \sqrt{3/5}$, $1$, and $1 + \sqrt{3/5}$ are fixed points. (b) Show that none of the three fixed points is 
locally convergent. 

Futher, find initial guesses for which FPI (c) cycles endlessly through numbers in the interval $(0,1)$, (d) the same as (c), but the interval is $(1, 2)$ (e) diverges to infinity. Cases (c) and (d) are examples of chaotic dynamics. In all three cases, FPI is unsuccessful.


```python
def f(x):
    return 1 - 6*x + 7.5*x**2 - 2.5*x**3

a = 0; b = 2

def g(x):
    return 1 - 5*x + 7.5*x**2 - 2.5*x**3


x0 = 0.5
x_star = 1 - np.sqrt(3./5)
w = interactive(Iteration_demo2, f=fixed(f), a=fixed(a), b=fixed(b), g=fixed(g), x0 = fixed(x0), iters=widgets.IntSlider(min=0,max=80,value=0))
display(w)
```


    interactive(children=(IntSlider(value=0, description='iters', max=80), Output()), _dom_classes=('widget-intera…



```python
x = fixed_iters(g, x0, 120)
plt.plot(np.arange(np.size(x)),x)
```


```python

```

db814818f4247e03c2c41a1b13734c2ddbd68a93
ipynb
Jupyter Notebook
Solution_of_Nonlinear.ipynb
xiaozhouli/numerical_analysis
68600ca56f8fdec6a2d22c65ec02ec2871546ea2
MIT
2021-11-06T02:11:33.000Z
2021-11-06T02:11:33.000Z
Solution_of_Nonlinear.ipynb
xiaozhouli/numerical_analysis
68600ca56f8fdec6a2d22c65ec02ec2871546ea2
MIT
Solution_of_Nonlinear.ipynb
xiaozhouli/numerical_analysis
68600ca56f8fdec6a2d22c65ec02ec2871546ea2
MIT
Qwen/Qwen-72B
1. YES
2. YES