Extended exercise

Compute the following integrals with the Midpoint rule, the Trapezoidal rule, and Simpson's rule: $$ \begin{eqnarray*} \int_{0}^{\pi}\sin x\, dx &=& 2,\\ \int_{-\infty}^{\infty} {1\over\sqrt{2\pi}} e^{-x^2}dx &=& 1,\\ \int_{0}^1 3x^2dx &=& 1,\\ \int_0^{\ln 11} e^xdx &=& 10,\\ \int_0^1 {3\over 2}\sqrt{x}dx &=& 1. \end{eqnarray*} $$ For each integral, write out a table of the numerical error for the three methods using a \( n \) function evaluations, where \( n \) varies as \( n=2^k+1 \), \( k=1,2, ..., 12 \).

A Python solution

In the extended problem, Solution 1 is obviously inferior because we need to apply, e.g., the Trapezoidal rule to five different integrand functions for 12 different \( n \) values. Then it only makes sense to implement the rule in a separate function that can be called 60 times.

Similarly, a mathematical function to be integrated is needed in three different rules, so it makes sense to isolate the mathematical formula for the integrand in a function in the language we are using.

We can briefly sketch a compact and smart Python code, in a single file, that solves the extended problem:

def f1(x):
    return sin(x)

def f2(x):
    return 1/sqrt(2)*exp(-x**2)

...

def f5(x):
    return 3/2.0*sqrt(x)

def Midpoint(f, a, b, n):
    ...

def Trapezoidal(f, a, b, n):
    ...

def Simpson(f, a, b, n):
    ...

problems = [(f1, 0, pi),  # list of (function, a, b)
            (f2, -5, 5),
            ...
            (f3, 0, 1)]

methods = (Midpoint, Trapezoidal, Simpson)
result = []
for method in methods:
    for func, a, b in problems:
        for k in range(1,13):
            n = 2**k + 1
            I = method(func, a, b, n)
            result.append((I, method.__name__, func.__name__, n))

# write out results, nicely formatted:
for I, method, integrand, n in result:
    print '%-20s, %-3s, n=%5d, I=%g' % (I, method, integrand, n)

Note that since everything in Python is an object that can be referred to by a variable, it is easy to make a list methods (list of Python functions), and a list problems where each element is a list of a function and its two integration limits. A nice feature is that the name of a function can be extracted as a string in the function object (name with double leading and trailing underscores).

To summarize, Solution 2 or 3 can readily be used to solve the extended problem, while Solution 1 is not worth much. In courses with many very simple exercises, solutions of type 1 will appear naturally. However, published solutions should employ approach 2 or 3 of the mentioned reasons, just to train students to think that this is a general mathematical method that I should make reusable through a function.

Along with the code above there should be a file test_integration_methods.py containing test functions for the various rules. The error formula for Simpson's rule contains \( f'''' \), so one can integrate a third-degree polynomial in a test and expect an error about the machine precision. The Midpoint rule integrates linear functions exactly. For testing of convergence rates, the Trapezoidal and Midpoint rules have errors behaving as \( n^{-2} \), while the error in Simpson's rule goes like \( n^{-4} \).

Solution 4: a Java OO program

Introductory courses in computer programming, given by a computer science department, often employ the Java language and emphasize object-oriented programming. Many computer scientists argue that it is better to start with Java than Python or (especially) Matlab. But how well is Java suited for introductory numerical programming?

Let us look at our first integration example, now to be solved in Java. Solution 1 is implemented as a simple main method in a class, with a code that follows closely the displayed Matlab code. However, students are in a Java course trained in splitting the code between classes and methods. Therefore, Solution 2 should be an obvious choice for a Java programmer. However, it is not possible to have stand-alone functions in Java, functions must be methods belonging to a class. This implies that one cannot transfer a function to another function as an argument. Instead one must apply the principles of object-oriented programming and implement the function argument as a reference to a superclass. To call the "function argument", one calls a method via the superclass reference. The code below provides the details of the implementation:

import java.lang.*;

interface Func {  // superclass for functions f(x)
    public double f (double x);  // default empty implementation
}

class f1 implements Func {
    public double f (double t)
    { return Math.exp(-Math.pow(t, 4)); }
}

class Trapezoidal {
    public static double integrate
           (Func f, double a, double b, int n)
    {
        double h = (b-a)/((double)n);
        double s = 0.5*(f.f(a) + f.f(b));
        int i;
        for (i = 1; i <= n-1; i++) {
            s = s + f.f(a+i*h);
        }
        return h*s;
    }
}

class MainProgram {
    public static void main (String argv[])
    {
        double a = -2;
        double b = 2;
        int n = 1000;
        double result = Trapezoidal.integrate(f, a, b, n);
        System.out.println(result);
    }
}

From a computer science point of view, this is a quite advanced solution since it relies on inheritance and true object-oriented programming. From a mathematical point of view, at least when compared to the Matlab and Python versions, the code looks unnecessarily complicated. Many introductory Java courses do not cover inheritance and true object-oriented programming, and without mastering these concepts, the students end up with Solution 1. On this background, one may argue that Java is not very suitable for implementing this type of numerical algorithms.