Fixed point iteration and cobweb plot

Question

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

• jonrsharpe: 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.

Answer 1

Your questions

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).

Additional comments

• 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):
 """An (arbitrary) quadratic mathematical function.
 It has fixpoints 0 and 0.5.
 >>> f(0)
 0
 >>> f(0.5)
 0.5
 """
 return 2*x*(1-x)
def g(x):
 """An (arbitrary) quadratic mathematical function.
 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())