# Fixed point iteration and cobweb plot

I'm using Python to find fixed points of a given function and then draw a cobweb plot to visualize it. Thanks to this question, I have the core of the code written and can accomplish the task, but I have a number of questions about improving the functionality of the code.

import numpy as np
import matplotlib.pyplot as plt

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

class fixpt(object):

def __init__(self, x):
self.x = x
self.f = f

def search(self, N = 1000, epsilon = 1.0E-100, store = False):
x = self.x
y = self.f(x)
n = 0
if store: values = [(x,y)]
while abs(y-x) >= epsilon and n < N:
x = f(x)
n += 1
y = f(x)
if store: values.append((x,y))
if store:
return y, values
else:
if n >= N:
return "No fixed point for given initial condition"
else:
return x, n, y
def plot(self):
res, points = self.search(store = True)
xaxis = np.linspace(0,1,1000)
plt.plot(xaxis, f(xaxis), 'b')
plt.plot(xaxis, xaxis, 'r')

for x, y in points:
plt.plot([x, x], [x, y], 'g')
plt.plot([x, y], [y, y], 'g')

plt.show()


What I'm wondering is how I can:

1. Avoid the need to redefine the function f(x) when analyzing other functions (by, perhaps asking for user input)
2. Avoid the need to change the range of concern (i.e. the range in np.linspace) in a similar manner
3. Improve the structure of this code in general (I'm a Python noob and get the feeling I've created a class for little to no reason
4. Most importantly: use this code iteratively to scan for fixed points for all x in a given range
• The answer to questions 1, 2 and 4 is to make the things you want to change into arguments to the functions/methods. And yes, there's no real reason for the class - search and plot could be standalone functions, they don't share any instance attributes. – jonrsharpe Jul 6 '14 at 19:50

1. Functions can be function arguments, too.
2. You can add a parameter for that.
3. See additional comments. I think it's fine that you use a class.
4. That is usually not in the range of codereview. That should be asked in StackOverflow (but I think I know what you want - just look at the code).

• Add a shebang. If you do, you can directly execute the script. See What does #!/usr/bin/python mean?
• Read PEP8 or at least use pep8online.com to make your coding style standard conform and hence easier to read. It's mostly about setting whitespaces correct.
• Use docstrings and eventually doctest. That makes your code again easier to read and you can easily verify if it works as you expect.
• If you don't use a variable, you should remove it at all. If that is not possible (e.g. in the case of unzipping values) you should name the variable _.

## The code

#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt

def f(x):
It has fixpoints 0 and 0.5.

>>> f(0)
0
>>> f(0.5)
0.5
"""
return 2*x*(1-x)

def g(x):
It has fixpoints 0 and (+/-) 1/sqrt(2)

>>> g(0)
0
>>> "%0.4f" % g(-0.7071)
'-0.7071'
"""
return 2*x*(1-x)*(1+x)

class fixpt(object):
"""Find fixed points of a given function f."""
def __init__(self, x, f, xmin=0, xmax=1, epsilon=1.0E-100):
"""Constructor
@param x float - where to start search
@param f function - function that should get evaluated
@param xmin float - Interval that gets examined
@param xmax float - Interval that gets examined
@param epsilon float - comparison when floats are considered to be
the same
"""
self.x = x
self.f = f
self.xmin = xmin
self.xmax = xmax
self.epsilon = epsilon

def search(self, N=1000, store=False):
"""Search fixed points.
@param N int - number of iterations
@param store - bool - if set, return ????
@return I did not understand that. It seems to depend on 'store'

>>> fixpointiter = fixpt(0, f)
>>> fixpointiter.search()
(0, 0, 0)
"""
# Set initial values
x = self.x
y = self.f(x)
n = 0
if store:
values = [(x, y)]  # Why?
while abs(y-x) >= self.epsilon and n < N:
x = f(x)
n += 1
y = f(x)
if store:
values.append((x, y))
if store:
return y, values
else:
if n >= N or x == -float('inf') or x == float('inf'):
# No fixed point for given initial condition
return None
else:
return x, n, y

def plot(self):
"""Draw a cobweb plot."""
_, points = self.search(store=True)
xaxis = np.linspace(self.xmin, self.xmax, 1000)
plt.plot(xaxis, f(xaxis), 'b')
plt.plot(xaxis, xaxis, 'r')

for x, y in points:
plt.plot([x, x], [x, y], 'g')
plt.plot([x, y], [y, y], 'g')

plt.show()

def scan(self):
"""Scan for all fix points in xmin / xmax."""
x_range = np.linspace(self.xmin, self.xmax, 1000)
fixpoints = []
for x in x_range:
self.x = x
ret = self.search()
if ret is not None:
fixpoints.append(ret[0])
# Get rid of duplicates
fixpoints = sorted(fixpoints)
final = []
for x in fixpoints:
# It's not the same as the last fixpoint
if len(final) == 0 or abs(final[-1]-x) >= self.epsilon:
final.append(x)
return final

if __name__ == '__main__':
# Automatic doctest
import doctest
doctest.testmod()
# Usage example
fixpointiter = fixpt(0, g, -2, 2)
print(fixpointiter.search())
print(fixpointiter.search(store=True))
print(fixpointiter.scan())