## Introduction

The Simpson’s rule is a numerical method used to approximate the definite integral of a function. It is a closed Newton-Cotes formula that uses quadratic polynomials to approximate the function on each subinterval of the partition. It is named after Thomas Simpson, an 18th-century mathematician who first described the method. The Simpson’s rule is widely used in numerical analysis, physics, and engineering to solve problems that involve integration.

## Implementation

The following implementation is in OCaml:

``````let simpson ?(n = 20) ?(eps = 1e-6) f a b =
let s_new = ref 0. in
let s_old = ref 0. in
let o_new = ref 0. in
let o_old = ref 0. in
(try
for i = 1 to n do
s_new := trapzd f a b i;
s_old := ((4. *. !s_new) -. !o_new) /. 3.;
if i > 5
then (
let d = abs_float (!s_old -. !o_old) in
let e = eps *. abs_float !o_old in
assert (not (d < e || (!s_old = 0. && !o_old = 0.)));
o_old := !s_old;
o_new := !s_new)
done
with
| _ -> ());
!s_new
``````

The function `simpson` takes in three arguments: a function `f` that returns a float, and two floats `a` and `b` representing the lower and upper bounds of the integral. It also has two optional parameters: `n` is the number of intervals used in the approximation, and `eps` is the maximum relative error allowed. The default values for `n` and `eps` are 20 and 1e-6, respectively.

## Step-by-step Explanation

1. The `simpson` function takes in three arguments: a function `f` that returns a float, and two floats `a` and `b` representing the lower and upper bounds of the integral. It also has two optional parameters: `n` is the number of intervals used in the approximation, and `eps` is the maximum relative error allowed. The default values for `n` and `eps` are 20 and 1e-6, respectively.

2. The function initializes four variables: `s_new`, `s_old`, `o_new`, and `o_old`, all of which are set to 0. `s_new` and `s_old` will be used to store the current and previous approximations of the integral, respectively, while `o_new` and `o_old` will be used to store the current and previous approximations of the error, respectively.

3. The function uses a `try`-`with` block to catch any exceptions that might be raised during the execution of the loop. This is done to ensure that the function always returns a value, even if an error occurs.

4. The function then enters a `for` loop that iterates from 1 to `n`. For each iteration of the loop, the function computes the current approximation of the integral over the interval [a,b] using the `trapzd` function, and stores it in the `s_new` variable.

5. The function then uses the current and previous approximations of the integral to compute a new approximation that is more accurate. This is done using the formula:

``````s_old = (4 * s_new - o_new) / 3
``````

where `s_old` is the new approximation, `s_new` is the current approximation, and `o_new` is the previous approximation.

1. The function then checks if the current iteration is greater than 5. If it is, the function checks if the relative error between the current and previous approximations is less than the maximum relative error allowed (`eps`). If it is, the function returns the current approximation of the integral (`s_new`). If it is not, the function updates the `o_old` and `o_new` variables to store the current and previous approximations of the error, respectively.

2. The function exits the loop and returns the current approximation of the integral (`s_new`).

3. If an exception is raised during the execution of the loop, the function catches it and returns the current approximation of the integral (`s_new`).

Overall, the `simpson` function uses the Simpson’s rule to approximate the integral of a function over an interval [a,b]. It does this by dividing the interval into `n` subintervals, approximating the function on each subinterval using a quadratic polynomial, and summing the integrals over each subinterval to obtain an approximation of the integral over the entire interval. The function uses the `trapzd` function to approximate the integral over each subinterval, and uses the current and previous approximations of the integral to compute a new approximation that is more accurate. The function also checks if the relative error between the current and previous approximations is less than the maximum relative error allowed, and returns the current approximation if it is.

## Complexity Analysis

The time complexity of the `simpson` function depends on the number of intervals used in the approximation (`n`) and the maximum relative error allowed (`eps`).

The `trapzd` function, which is used to approximate the integral over each subinterval, has a time complexity of O(n), where n is the number of intervals used in the approximation. This is because the function uses a `for` loop that iterates `n` times to compute the sum of the areas of the trapezoids.

The `simpson` function uses the `trapzd` function `n` times to approximate the integral over each subinterval, and then uses the current and previous approximations of the integral to compute a new approximation that is more accurate. The function also checks if the relative error between the current and previous approximations is less than the maximum relative error allowed, and returns the current approximation if it is.

The time complexity of the `simpson` function can be approximated as follows:

• The `for` loop that iterates from 1 to `n` has a time complexity of O(n).
• The `trapzd` function is called `n` times, so its time complexity is O(n^2).
• The time complexity of the computations inside the loop is O(1).
• The `if` statement that checks if the relative error is less than the maximum relative error allowed has a time complexity of O(1).
• The `try`-`with` block has a time complexity of O(1).
• The function returns the current approximation of the integral, which has a time complexity of O(1).

Therefore, the overall time complexity of the `simpson` function is O(n^2). This means that the function takes longer to run as the number of intervals used in the approximation (`n`) increases. The maximum relative error allowed (`eps`) does not affect the time complexity of the function, but it does affect the accuracy of the approximation. A smaller value of `eps` will result in a more accurate approximation, but will also require more iterations of the `for` loop and hence increase the running time of the function.