direction_field.py

import numpy as np
from matplotlib import pyplot


class Viewport:
    def __init__(self, rng, *rest):
        if len(rest) == 0:
            self.array = [-rng, rng, -rng, rng]
        elif len(rest) == 1:
            self.array = [rng, rest[0]] * 2
        elif len(rest) == 3:
            self.array = [rng] + list(rest)
        else:
            raise Exception('Viewport() only accepts 1, 2, or 4 arguments')


class DiffEqPlot:
    def __init__(self, axis):
        self.plot = axis
        self.plot.set_title(r'Direction Fields', fontsize=20)
        self.plot.set_xlabel(r't', fontsize=15)
        self.plot.set_ylabel(r'f(t)', fontsize=15)
    def display(self, view, dif_eq, step=0.4, color='r'):
        self.plot.axis(view.array)
        dt = step

        for t in np.arange(view.array[0], view.array[1], step):
            for f in np.arange(view.array[2], view.array[3], step):
                self.plot.plot([t], [f], 'ko', markersize=2.0)
                if f == 0:
                    f = 0.00000001
                if t == 0:
                    t = 0.00000001
                df = dif_eq(f, t) * dt
                mag = np.sqrt(df**2 + dt**2) / dt
                self.plot.arrow(t, f, dt / mag, df / mag, head_width=0.04,
                head_length=0.08, fc='b', ec=color)


if __name__ == '__main__':
    #pyplot.rc('text', usetex=True)
    #pyplot.rc('font', family='serif')

    ax = pyplot.subplot()

    plot = DiffEqPlot(ax)
    plot.display(Viewport(10), lambda f, t: (-f*(1+t))/(t+(2/np.e**t)))
    pyplot.show()

euler.py

import numpy as np

def estimate_diff_eq(de, interval, step_size, initial_value):
    f = initial_value
    n_steps = ((interval[1] - interval[0]) / step_size) + 1

    for t in np.linspace(interval[0], interval[1], n_steps):
        df = de(f, t)
        yield (t, f, df)
        f += df * step_size


if __name__ == '__main__':
    def diff_eq(f, t):
        return f**2/(t-2)

    table = list(estimate_diff_eq(diff_eq, (0, 3), 0.75, 4))

    for e in table:
        print('{:.2f}\t{:.2f}\t{:.2f}'.format(e[0], e[1], e[2]))