Integrata: Numerical Integration Demo (Julia & Fortran)

Integrata: Numerical Integration Demo (Julia & Fortran)#

Welcome to the Integrata documentation and demo notebook! This notebook demonstrates the concepts behind Integrata using the Julia language, with code and explanations for Simpson’s Rule and the Trapezoidal Rule. We’ll also compare results to Python/NumPy/Scipy and show how to use Integrata from the command line.

1. Import Required Libraries#

We’ll use Julia’s built-in math and plotting libraries for demonstration.

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
print("Libraries loaded. Ready for integration demos!")

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 2
      1 import numpy as np
----> 2 import matplotlib.pyplot as plt
      4 np.random.seed(42)
      5 print("Libraries loaded. Ready for integration demos!")

ModuleNotFoundError: No module named 'matplotlib'

2. Simpson’s Rule: Theory & Example#

Simpson’s Rule is a numerical method for approximating definite integrals. It works best for smooth functions and divides the interval into an even number of subintervals.

We’ll demonstrate Simpson’s Rule by integrating \(f(x) = \sin(x)\) over \([0, \pi]\), where the exact result is 2.

import numpy as np
import matplotlib.pyplot as plt

f = np.sin
a, b = 0.0, np.pi

xs = np.linspace(a, b, 100)
ys = f(xs)

plt.plot(xs, ys, label="sin(x)")
plt.xlabel("x")
plt.ylabel("f(x)")
plt.title("f(x) = sin(x) on [0, π]")
plt.legend(loc="upper right")
plt.show()

def simpsons_rule(f, a, b, n):
    if n % 2 != 0:
        raise ValueError("n must be even for Simpson's Rule")
    h = (b - a) / n
    s = f(a) + f(b)
    for i in range(1, n, 2):
        s += 4 * f(a + i * h)
    for i in range(2, n-1, 2):
        s += 2 * f(a + i * h)
    return s * h / 3

n = 10  # Must be even
approx = simpsons_rule(f, a, b, n)
print(f"Simpson's Rule approximation (n={n}): {approx:.8f}")
print(f"Exact value: {2.0:.8f}")
print(f"Absolute error: {abs(approx - 2.0):.2e}")

3. Trapezoidal Rule: Theory & Example#

The Trapezoidal Rule approximates the area under a curve by dividing it into trapezoids. It’s simpler but less accurate than Simpson’s Rule for smooth functions.

Let’s apply it to \(f(x) = \sin(x)\) over \([0, \pi]\) for comparison.

def trapezoidal_rule(f, a, b, n):
    h = (b - a) / n
    s = 0.5 * (f(a) + f(b))
    for i in range(1, n):
        s += f(a + i * h)
    return s * h

n = 10
trap_approx = trapezoidal_rule(f, a, b, n)
print(f"Trapezoidal Rule approximation (n={n}): {trap_approx:.8f}")
print(f"Exact value: {2.0:.8f}")
print(f"Absolute error: {abs(trap_approx - 2.0):.2e}")

4. Results Comparison & Using Integrata#

As shown above, Simpson’s Rule is generally more accurate than the Trapezoidal Rule for smooth functions like \(\sin(x)\). Both methods are implemented in Integrata (Fortran), which you can use from the command line:

make run

Follow the menu prompts to select the integration method and input your function, interval, and number of subintervals.

For more details, see the README and the Fortran source code in src/.

5. Further Reading & Contributing#

Want to contribute? See CONTRIBUTING.md and open an issue or pull request!

This notebook is part of the Integrata documentation and can be published on GitHub Pages using Jupyter Book or similar tools.

Contents