<a href="https://colab.research.google.com/github/stephenbeckr/numerical-analysis-class/blob/master/Demos/Ch2_FixedPointPlots.ipynb" target="_parent"></a>

# Fixed Point iteration (and cobweb plots)

See [cobweb plot on wikipedia](https://en.wikipedia.org/wiki/Cobweb_plot), sometimes called a spider diagram

You can plot these via Python, but to make them nice and interactive is a bit of work (`matplotlib.animation.FuncAnimation` can do it, but you have to be quite proficient with Python and matplotlib).

For this demo, I recommend that you use [this geogebra online cobweb plot](https://www.geogebra.org/m/uvsfvNDt) which is quite simple and runs on most modern web browsers.

We'll look at 5 different fixed-point reformulations of the same root-finding problem
$$f(x) = x^3 + 4x^2 - 10, \text{ find } f(x)=0$$

where we now solve the fixed point equation $$g(x)=x$$ for:
1. $$g_1(x) = x - x^3 - 4x^2+10$$
2. $$g_2(x) = \sqrt{\frac{10}{x}-4x}$$
3. $$g_3(x) = \frac{1}{2}\sqrt{10-x^3}$$
4. $$g_4(x) = \sqrt{\frac{10}{4+x}}$$
5. $$g_5(x) = x - \frac{x^3+4x^2-10}{3x^2+8x}$$

These are in the book in Table 2.2.

The goal of this demo is:
1. Learn how to use and interpret a cobweb plot
2. Predict which of the fixed point formulations will converge (hint: plot their deriviative)
3. Check your predictions with the cobweb plot


```
import math
from sympy import init_printing
init_printing()

#from sympy import symbols
from sympy import *
from sympy.plotting import plot

x = symbols('x')
f, fprime = symbols('f fprime',cls=Function)

f = x**3 + 4*x**2 - 10
plot( f, (x,0,3), ylim = (-5,5) );
```


```
from scipy.optimize import root_scalar
output = root_scalar( lambda x :x**3 + 4*x**2 - 10, x0 = 1.5, x1 = 1.4)
output
```




          converged: True
               flag: 'converged'
     function_calls: 6
         iterations: 5
               root: 1.3652300134140967




```
# The fix point iteration is very simple.  Not having any error checking or convergence checks, 
# in order to focus on the algorihtm itself, the following is an implementation
def fixedPointIter( f, x0, maxIter ):
  p = x0 # if a numpy array, then do x0.copy()
  for i in range(maxIter):
    p = f(p)
    print("Iter {:3d}, guess is {:12.9f}".format(i,p))
  return p

p = fixedPointIter( lambda x : .5*math.sqrt(10-x**3), 1.5, 10 )
p
```

## Formulation 1
$$g_1(x) = x - x^3 - 4x^2+10$$


```
g1 = x - x**3 - 4*x**2+10
g1
```


```
g1Prime = diff(g1)
g1Prime
```


```
plt1 = plot( g1, (x,0,3), ylim = (-5,5), show=False )
plt2 = plot( g1Prime, (x,0,3), show=False , line_color='red')
plt3 = plot( x, (x,0,3), show=False , line_color='green')
plt1.append( plt2[0] )
plt1.append( plt3[0] )
#ax.axis('equal') # equal aspect ratio, so 45 degree line looks like 45 degrees
plt1.show()
```


```

```

## Formulation 2
$$g_2(x) = \sqrt{\frac{10}{x}-4x}$$


```
g2 = sqrt( 10/x-4*x )
g2
```


```
g2Prime = diff(g2)
g2Prime
```


```
plt1 = plot( g2, (x,1,2), ylim = (-10,3), show=False )
plt2 = plot( g2Prime, (x,1,2), show=False , line_color='red')
plt3 = plot( x, (x,1,2), show=False , line_color='green')
plt1.append( plt2[0] )
plt1.append( plt3[0] )
#ax.axis('equal') # equal aspect ratio, so 45 degree line looks like 45 degrees
plt1.show()
```

## Formulation 3
$$g_3(x) = \frac{1}{2}\sqrt{10-x^3}$$


```
g3 = 1/2*sqrt(10-x**3)
g3
```


```
g3Prime = diff(g3)
g3Prime
```


```
plt1 = plot( g3, (x,1,3), ylim = (-5,5), show=False )
plt2 = plot( g3Prime, (x,1,3), show=False , line_color='red')
plt3 = plot( x, (x,1,3), show=False , line_color='green')
plt1.append( plt2[0] )
plt1.append( plt3[0] )
plt1.show()
```

## Formulation 4
$$g_4(x) = \sqrt{\frac{10}{4+x}}$$


```
g4 = sqrt(10/(4+x))
g4
```


```
g4Prime = diff(g4)
g4Prime
```


```
plt1 = plot( g4, (x,1,3), ylim = (-5,5), show=False )
plt2 = plot( g4Prime, (x,1,3), show=False , line_color='red')
plt3 = plot( x, (x,1,3), show=False , line_color='green')
plt1.append( plt2[0] )
plt1.append( plt3[0] )
plt1.show()
```

## Formulation 5
$$g_5(x) = x - \frac{x^3+4x^2-10}{3x^2+8x}$$


```
g5 = x - (x**3+4*x**2-10)/(3*x**2+8*x)
g5
```


```
g5Prime = diff(g4)
g5Prime
```


```
plt1 = plot( g5, (x,1,3), ylim = (-5,5), show=False )
plt2 = plot( g5Prime, (x,1,3), show=False , line_color='red')
plt3 = plot( x, (x,1,3), show=False , line_color='green')
plt1.append( plt2[0] )
plt1.append( plt3[0] )
plt1.show()
```


```

```

9bd9a35148d505d8da3cc9d0d4021783d734d5a5
ipynb
Jupyter Notebook
Demos/Ch2_FixedPointPlots.ipynb
skhadem/numerical-analysis-class
a022fc2e73254800a1e193f94280446223ef01ae
BSD-3-Clause
2020-08-25T19:11:46.000Z
2021-11-12T21:49:44.000Z
Demos/Ch2_FixedPointPlots.ipynb
sabrinazhengliu/numerical-analysis-class
1585b68d3f3c375259507bf18f03ba0332851edb
BSD-3-Clause
2020-09-01T21:44:12.000Z
2020-09-01T21:44:12.000Z
Demos/Ch2_FixedPointPlots.ipynb
sabrinazhengliu/numerical-analysis-class
1585b68d3f3c375259507bf18f03ba0332851edb
BSD-3-Clause
2020-08-25T21:25:17.000Z
2021-03-03T04:14:26.000Z
Qwen/Qwen-72B
1. YES
2. YES