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 \).
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} \).
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.