recursive_bessel.py

from scipy import special
from scipy import misc
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages


def relative_error(accepted, approximation):
    # Helper function for finding relative error.
    return abs((accepted - approximation) / float(accepted))


def err_of_next_bessel_iter(order, x, tolerance):
    '''
    For each call of err_of_next_bessel_iter(), this function takes another
    recursive step and returns the relative error of the current approximation
    to the modified bessel function specified by the arguments order and x. The
    iteration stops once the relative error is less than the specified
    tolerance.
    Since the modified Bessel function is a sum over m, at each iteration
    this function keeps track of the previous term in order to calculate the
    current term, and adds the current term to the running sum.

    order: int          the order of the modified bessel function
    x: float            the position at which to evaluate the modified bessel
                        function
    tolerance: float    the tolerance for which the approximation is "good
                        enough"
    '''
    # Set up the variables for the m = 0 case.
    accepted_value = special.iv(order, x)
    current_term = (1.0 / misc.factorial(order)) * (x / 2.0) ** order
    current_sum = current_term
    current_error = relative_error(accepted_value, current_sum)
    m = 0

    # n digits of accuracy requires a tolerance of 10 ** (-n)
    while current_error > tolerance:
        # The next term in the summation is the current term with
        # modifications to the two factorials and the order of (x/2).
        next_term = current_term * \
                    ((x / 2.0) ** 2) / ((m + 1) * (m + order + 1))
        current_term = next_term
        current_sum += current_term
        m += 1

        current_error = relative_error(accepted_value, current_sum)
        yield current_error

if __name__ == '__main__':
    xs = range(1, 11)
    data = [[] for _ in xrange(10)]
    for i, x in enumerate(xs):
        # The list() function stores the output of all the iterations of
        # err_of_next_bessel_iter() in a list.
        data[i] = np.array(list(err_of_next_bessel_iter(2, x, 10 ** (-8))))
    plots = np.array([(np.arange(1, len(d) + 1), d) for d in data])
    plots = plots.flatten()

    _, ax = plt.subplots()
    plt.grid()
    plt.title("Precision of I_2(x) Calculation v. Number of Steps", size=18)
    plt.xlabel("Number of Back Recursion Steps", size=15)
    x_right_lim = max(map(len, plots))
    plt.xlim(1, x_right_lim)
    ax.xaxis.set_ticks(range(1, x_right_lim + 1))
    plt.ylabel("Relative Error of Calculation", size=15)
    ax.set_yscale('log')
    ax.plot(*plots)
    plt.legend(['I_2({})'.format(x) for x in xs], fontsize=13, ncol=2)
    plt.axhline(y=10 ** (-4))
    plt.axhline(y=10 ** (-8))
    plt.tight_layout()
    output = PdfPages('problem_3_graph.pdf')
    output.savefig()
    output.close()