Approximate methods for differential equations and integration

These examples are adapted from Sage wiki:

In [1]:
x,y = var('x,y')

Approximate integration

The following two cells implement three integration methods: midpoint, trapezoidal, Simpson's rules. In all the three methods the interval $[-1,1]$ is partitioned into $n$ subintervals. The corresponding sums consist of n (midpoint) and n+1 terms (trapezoidal, Simpson's rule), respectively.

In [2]:
import scipy
import numpy
from scipy.special.orthogonal import p_roots, t_roots, u_roots
from scipy.integrate import quad, trapz, simps
from sage.ext.fast_eval import fast_float
from numpy import linspace

methods = {
    'Midpoint': {'w': 1, 'xmin': -1, 'xmax': 1, 
        'func': lambda n: (linspace(-1r+1.0r/n,1r-1.0r/n,n),
                'nterms': lambda n: n},
     'Trapezoid': {'w': 1, 'xmin': -1, 'xmax': 1, 
        'func': lambda n: (linspace(-1r,1r,n+1),
                  'nterms': lambda n: n+1},
     'Simpson': {'w': 1, 'xmin': -1, 'xmax': 1, 
        'func': lambda n: (linspace(-1r,1r,n+1), 
                'nterms': lambda n: n}}
def box(center, height, width2,**kwds):
    return polygon([(center-width2,0),
In [3]:
def weights(n=slider(2,30,2,default=10),
    ff = fast_float(f,'x')
    method = methods[show_method]
    xcoords,w = (method['func'])(int(n))
    xmin = method['xmin']
    xmax = method['xmax']
    plot_min = max(xmin, -10)
    plot_max = min(xmax, 10)

    coords = zip(xcoords,w)
    max_weight = max(w)
    coords_scaled = zip(xcoords,w/max_weight)

    f_graph = plot(f,plot_min,plot_max)
#     boxes = sum(box(x,ff(x),(xmax-xmin)/n/2,rgbcolor=(0.5,0.5,0.5),alpha=0.3) for x,w in coords)
    stems = sum(line([(x,0),(x,ff(x))],rgbcolor=(1-y,1-y,1-y),
        thickness=2,markersize=6,alpha=y) for x,y in coords_scaled)
    points = sum([point([(x,0),
        (x,ff(x))],rgbcolor='black',pointsize=30) for x,_ in coords])
    graph = stems+points+f_graph#+boxes
    approximation = sum([w*ff(x) for x,w in zip(xcoords,w)])
    integral,integral_error = scipy.integrate.quad(ff, xmin, xmax)
    x_val = linspace(-1, 1,n)
    y_val = [*map(ff,x_val)]
    trapezoid = integral-trapz(y_val, x_val)
    simpson = integral-simps(y_val, x_val)
    pretty_print(html(r"$$\sum_{i=1}^{i=%s} w_i f(x_i^*) \Delta x = %s\approx \int_{-1}^{1}%s \,dx = %s $$"%
                      ( (method['nterms'])(int(n)), approximation,  latex(f), integral)))
    error_data = [integral-approximation, trapezoid, simpson, integral_error]
    print("Errors for integration methods:\n")
    print("Our implementation:\n%s\t %s\n\nBuilt-in:\nquad\t\t %s\ntrapz\t\t %s\nsimps\t\t %s" % (

In the formula above, $w_i$ are the weights on values of $f$, which show up in the integration rules. For example, we have seen that the n-th sum in the midpoint rule is $$ M_n = \Delta x\left(f(\bar x_1) + f(\bar x_2)+\ldots+ f(\bar x_n) \right). $$ We interpret this expression as the sum $$ \sum_{i=1}^n w_i f(\bar x_i)\,\Delta x, $$ where all the $w_i$ are equal to 1. On the other hand, for the trapezoidal rule, the n-th sum is $$ M_n = \frac{\Delta x}{2} \left(f(x_0) + 2f(\bar x_1)+\ldots+ 2f(x_{n-1}) + f(x_{n}) \right), $$ so in this sum $$ \sum_{i=0}^n w_i f(x_i)\,\Delta x, $$ $w_i$ are equal to $$ \frac 12, \ 1, \ 1, \ 1, \ldots, 1,\ 1, \ \frac12. $$ (Just multiply the 1/2 in front with the coefficients at $f(x_i)$.)

Simple case of Euler's method

The following cell implements Euler's method and compares the approximate solution obtained with it (blue) to the exact solution (red).

In [4]:
def euler_method(
    y_exact_in = input_box('exp(x)', type = str, label = 'Exact solution = '),
    y_prime_in = input_box('exp(x)', type = str, label = "y' = "),
    start = input_box(0.0, label = 'x starting value: '),
    startval = input_box(1.0, label = 'y starting value: '),
    stop = input_box(3.0, label = 'x stopping value: '),
    nsteps = slider([2^m for m in range(0,10)], default = 10, label = 'Number of steps: '),
    y_exact = lambda x: eval(y_exact_in)
#     dy/dx = F(x,y)
    y_prime = lambda x,y: eval(y_prime_in)
    stepsize = float((stop-start)/nsteps)
    sol = [startval]+[0.0]*nsteps
    xvals = [start]+[0.0]*nsteps
    for step in range(nsteps):
        sol[step+1] = sol[step] + stepsize*y_prime(xvals[step],sol[step])
        xvals[step+1] = xvals[step] + stepsize
    sol_max = max(sol + [find_local_maximum(y_exact,start,stop)[0]])
    sol_min = min(sol + [find_local_minimum(y_exact,start,stop)[0]])
    show(plot(y_exact(x),start,stop,rgbcolor=(1,0,0))+line([[xvals[index],sol[index]] for index in range(len(sol))]),xmin=start,xmax = stop, ymax = sol_max, ymin = sol_min)

Direction field compared to Euler's method

In the next two cells, we plot the vector field for the vector with coordinates $(f,g)$, where both $f$ and $g$ depend on $(x,y)$. The vector field puts at every point $(x,y)$ of the plain vector $$ (f(x,y), g(x,y).$$ In the setting from the lecture, where we had $$ \frac{dy}{dx} = F(x,y), $$ this means that $$ F(x,y) = \frac{g(x,y)}{f(x,y)}. $$ That is, the slope of a vector is given by the ratio of second coordinate to the first one.

The blue curve connects the points $(x_k, y_k)$ for $k\geq 0$, obtained by the iterations of the Euler's method. By making the step size ($\Delta x$ in the lecture) sufficiently small, we get a good approximation to the solution of the equation, corresponding to the vector field.

In [5]:
from sage.ext.fast_eval import fast_float
In [6]:
def _(f = input_box(default=y), g=input_box(default=-x*y+x^3-x),
      xmin=input_box(default=-1), xmax=input_box(default=1),
      ymin=input_box(default=-1), ymax=input_box(default=1),
      start_x=input_box(default=0.5), start_y=input_box(default=0.5),
      step_size= slider([2^m for m in range(-10,1)], default =0.125, label = r'$\Delta x$ '),
      steps=slider([2^m for m in range(0,11)], default =128, label = 'Number of steps: ') ):
    ff = fast_float(f, 'x', 'y')
    gg = fast_float(g, 'x', 'y')
    steps = int(steps)

    points = [ (start_x, start_y) ]
    for i in range(steps):
        xx, yy = points[-1]
            points.append( (xx + step_size * ff(xx,yy), yy + step_size * gg(xx,yy)) )
        except (ValueError, ArithmeticError, TypeError):

    starting_point = point(points[0], pointsize=50)
    solution = line(points)
    vector_field = plot_vector_field( (f,g), (x,xmin,xmax), (y,ymin,ymax) )

    result = vector_field + starting_point + solution

            r"$\displaystyle\frac{dy}{dx} = \frac gf = \frac{%s}{%s}$  <br> " % (latex(g),latex(y))
In [ ]:
In [ ]: