# Nonlinear Equations

We want to find a root of the nonlinear function $f$ using different methods.

1. Bisection method
2. Newton method
3. Chord method
4. Secant method
5. Fixed point iterations






```python
%matplotlib inline
from numpy import *
from matplotlib.pyplot import *
import sympy as sym

```


```python
t = sym.symbols('t')

f_sym = t/8. * (63.*t**4 - 70.*t**2. +15.) # Legendre polynomial of order 5

f_prime_sym = sym.diff(f_sym,t)

f = sym.lambdify(t, f_sym, 'numpy')
f_prime = sym.lambdify(t,f_prime_sym, 'numpy')

phi = lambda x : 63./70.*x**3 + 15./(70.*x)
#phi = lambda x : 70.0/15.0*x**3 - 63.0/15.0*x**5
#phi = lambda x : sqrt((63.*x**4 + 15.0)/70.)

# Let's plot
n = 1025

x = linspace(-1,1,n)
c = zeros_like(x)

_ = plot(x,f(x))
_ = plot(x,c)
_ = grid()

```


```python
# Initial data for the variuos algorithms

# interval in which we seek the solution 
a = 0.7
b = 1.

# initial points
x0 = (a+b)/2.0
x00 = b

```


```python
# stopping criteria
eps = 1e-10
n_max = 1000
```

## Bisection method

$$
x^k = \frac{a^k+b^k}{2}
$$
```
                   if (f(a_k) * f(x_k)) < 0:
                      b_k1 = x_k
                      a_k1 = a_k
                   else:
                      a_k1 = x_k
                      b_k1 = b_k
```


```python
def bisect(f,a,b,eps,n_max):
    assert f(a)*f(b)<0   # Check if initial values are valid
    a_new = a
    b_new = b
    x = mean([a,b])
    err = eps + 1.        # Make sure err > eps holds for first iteration
    errors = [err]
    it = 0
    while (err > eps and it < n_max):
        if ( f(a_new) * f(x) < 0 ):
            # root in (a_new,x)
            b_new = x
        else:
            # root in (x,b_new)
            a_new = x
        
        x_new = mean([a_new,b_new])
        
        #err = 0.5 *(b_new -a_new)
        err = abs(f(x_new))   # Use residual as error estimate
        #err = abs(x-x_new)
        
        errors.append(err)
        x = x_new
        it += 1
    
    semilogy(errors)
    print (it)
    print (x)
    print (err)
    return errors
              
errors_bisect = bisect(f,a,b,eps,n_max)


             
        
```


```python
# is the number of iterations coherent with the theoretical estimation?
```

In order to find out other methods for solving non-linear equations, let's compute the Taylor's series of $f(x^k)$ up to the first order 

$$
f(x^k) \simeq f(x^k) + (x-x^k)f^{\prime}(x^k)
$$
which suggests the following iterative scheme
$$
x^{k+1} = x^k - \frac{f(x^k)}{f^{\prime}(x^k)}
$$

The following methods are obtained applying the above scheme where

$$
f^{\prime}(x^k) \approx q^k
$$

## Newton's method
$$
q^k = f^{\prime}(x^k)
$$

$$
x^{k+1} = x^k - \frac{f(x^k)}{q^k}
$$


```python
def newton(f,f_prime,x0,eps,n_max):
    x = x0
    err = abs(f(x0))
    errors = [err]
    it = 0
    while (err>eps and it<n_max):
        qk = f_prime(x)
        if abs(qk) < 1e-12:
            raise RuntimeError("f_prime(x) is close to zero")
        x_new = x - f(x) / qk
        
        #err = abs(f(x_new))
        err = abs(x-x_new)
        errors.append(err)
        it += 1
        
        x = x_new
        
    print (it)
    print (err)
    semilogy(errors)

    return errors

%time errors_newton = newton(f,f_prime,1.,eps,n_max)
#x_new, errors = newton(f,f_prime,1.,eps,n_max)
#print(x_new)

```

## Chord method

$$
q^k \equiv q = \frac{f(b)-f(a)}{b-a}
$$

$$
x^{k+1} = x^k - \frac{f(x^k)}{q}
$$


```python
def chord(f,a,b,x0,eps,n_max):
    x = x0
    err = abs(f(x0))
    it = 0
    q = (f(b)-f(a))/(b-a)  # use estimation for q instead of real derivative
    
    errors=[err]
    while (err>eps and it<n_max):
        x_new = x - f(x) / q
        
        # err = abs(f(x_new))
        err = abs(x_new-x)
        errors.append(err)
        it += 1
        
        x = x_new
    print (it)
    print (err)
    semilogy(errors)

    return errors

errors_chord = chord (f,a,b,x0,eps,n_max)
```

## Secant method

$$
q^k = \frac{f(x^k)-f(x^{k-1})}{x^k - x^{k-1}}
$$

$$
x^{k+1} = x^k - \frac{f(x^k)}{q^k}
$$

Note that this algorithm requirs **two** initial points


```python
def secant(f,x0,x00,eps,n_max):
    x = x0
    x_new = x00
    err = abs(f(x0))
    errors = [err]
    it = 0
    while (err>eps and it<n_max):
        qk = (f(x_new)-f(x))/(x_new-x)
        if abs(qk) < 1e-12:
            raise RuntimeError("q_k is close to zero")
        x = x_new  # store actual value before overwriting
        x_new = x_new - f(x_new) / qk
        
#        err = abs(f(x_new))
        err = abs(x_new-x)
        errors.append(err)
        it += 1
    
    print (it)
    print (err)
    semilogy(errors)  
    
    return errors
    

errors_secant = secant(f,x0,x00,eps,n_max)
```

## Fixed point iterations

$$
f(x)=0 \to x-\phi(x)=0
$$

$$
x^{k+1} = \phi(x^k)
$$


```python
def fixed_point(phi,x0,eps,n_max):
    x = x0
    
    err = abs(phi(x0)-x0)
    errors=[err]
    it = 0
    while (err>eps and it<n_max):
        x_new = phi(x)
        
        err = abs(phi(x)-x)  # residual
        err = abs(x_new-x)
        errors.append(err)
        
        x = x_new
        
    semilogy(errors)
    
    return errors

errors_fixed = fixed_point(phi,0.3,eps,n_max)
        
```

## Comparison


```python
# plot the error convergence for the methods
loglog(errors_bisect, label='bisect')
loglog(errors_chord, label='chord')
loglog(errors_secant, label='secant')
loglog(errors_newton, label ='newton')
loglog(errors_fixed, label ='fixed')
_ = legend()
```


```python
# Let's compare the scipy implmentation of Newton's method with our..
```


```python
import scipy.optimize as opt
%time opt.newton(f, 1.0, f_prime, tol = eps)
```

    CPU times: user 402 µs, sys: 6 µs, total: 408 µs
    Wall time: 417 µs





    0.906179845938664




```python

```

1addeeef3bec12bf33e4ced66aa46ed339a5d124
ipynb
Jupyter Notebook
notes/02_non_linear_equations.ipynb
luschu/P1.4_seed
1c44f6f54f0544faac97cc73ddd447f652a1ef71
CC-BY-4.0
notes/02_non_linear_equations.ipynb
luschu/P1.4_seed
1c44f6f54f0544faac97cc73ddd447f652a1ef71
CC-BY-4.0
notes/02_non_linear_equations.ipynb
luschu/P1.4_seed
1c44f6f54f0544faac97cc73ddd447f652a1ef71
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES