HiveBrain v1.2.0
Get Started
← Back to all entries
debugpythonMinor

Fixed point iteration and cobweb plot

Submitted by: @import:stackexchange-codereview··
0
Viewed 0 times
pointplotanditerationcobwebfixed

Problem

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:
                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:

  • Avoid the need to redefine the function f(x) when analyzing other functions (by, perhaps asking for user input)



  • Avoid the need to change the range of concern (i.e. the range in np.linspace) in a similar manner



  • 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



  • Most importantly: use this code iteratively to scan for fixed points for all x in a given range

Solution

Your questions

  • Functions can be function arguments, too.



  • You can add a parameter for that.



  • See additional comments. I think it's fine that you use a class.



  • 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 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())

Code Snippets

#!/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__

Context

StackExchange Code Review Q#56281, answer score: 4

Revisions (0)

No revisions yet.