Implementation of cumulative simpson method

Question

I am working on a program which uses cumulative integration methods to solve differential equations, and I want to confirm that my implementation of the Simpson method is correct, as I could not find any implementations of the Simpson method in a cumulative fashion online (except in a suspicious looking matlab file).

Function

def integrate(tt_step, new_val, integral_val):
    if step == 0:
        return tt_step*new_val
    elif step == max_steps:
        return (tt_step*new_val)+integral_val
    else:
        if step % 2 == 0:
            return (2*tt_step*new_val) + integral_val
        else:
            return (4*tt_step*new_val) + integral_val

(Written in Python 3.)

Description

Basically, the function takes as inputs tt_step which is the time step divided by 3, which I considered equivalent to $\frac{\Delta x}{3}$ (my function is being used to find the movement of a charged ion in 2-d space, so all integration is wrt time), new_val, which is the next function value, which I considered equivalent to $f(x_n)$, and integral_val, which is the integral from the previous time - that is, the integral up through $f(x_{n-1})$ where $n$ is the time.

The function first checks if step, or the time, is equal to zero, in which case it multiplies the time step divided by three, tt_step, by new_val; this basically returns zero. It then checks if the time is equal to the time the user wants to the program to end, max_steps, in which case it returns tt_step * new_val + integral_val. The reason these two cases are pulled out from the rest of the code is because these two cases are the only ones where new_val does not need to be multiplied by some coefficient ($2$ or $4$).

The rest of the code then checks if $t\mod 2 \equiv 0$ (where $t$ is time is step) and if so multiplies new_val by $2$ as well, and if not, then multiplies by $4$ as well. My reasoning for this is as follows: step starts at $0$. Then:

• At step $0$, $f(x_0)\rightarrow 1$
• At step $1$, $f(x_1)\rightarrow 4$
• At step $2$, $f(x_2)\rightarrow 2$
• At step $3$, $f(x_3)\rightarrow 4$
• At step $4$, $f(x_4)\rightarrow 2$
• $...$

Finally, tt_step is multiplied only by new_val in each step, not integral_val, because it has already been multiplied by integral_val.

Calling the function

To call the function (let us say we are trying to find the x position of the particle), we create a variable x_position_integral which we set equal to the initial position of the particle. We then have a loop, which starts at step = 0 and ends at step = max_steps. In this loop, we have a variable equal to the x velocity at that time, which updates with each loop. Thus we call the function each loop as x_position_integral = integrate(tt_step, x_velocity, x_position_integral and append each instance of x_position_integral to a list. The list is then graphed at the end of the program.

Question

Is this an accurate implementation of a cumulative version of the Simpson method? If there are any details that are unclear about the program, please let me know. To be clear, this is not so much about the programming, this is about whether or not this is mathematically correct.

Answer 1

This is rather wrong. (Sorry, me.)

As proof, consider this dandy little spreadsheet, with your formula applied to the function $x^2$ and a normal calculation of the integral. The average error over fifty steps with your function is 429.17. Eesh. (Ignore the "trap" column and the error column just to the right of it, as well as the "mod" column.)

As a general summary of what's wrong - you're summing over the course of the whole program, instead of doing that sum each step. The real code is more along the lines of:

def simp_integrate_correctly(t_step, new_val, integral_val):
    b_a = new_val[-1] - new_val[-2]
    return t_step/12 * (new_val[-2] + 4*(new_val[-2] + b_a/4) + 2*(new_val[-2] + 2*b_a/4) + 4*(new_val[-2] + 3*b_a/4) + new_val[-1]) + integral_val

And behold - again in the dandy spreadsheet, with this new formula applied to the function $x^2$, the average error is -0.000000007119216775 over 500 steps.

Much better.

Answer 2

It looks like you're on the right track, but hard to make a full assessment without seeing a bit more code. For example, where is the loop over step?

Here's a link to the Wikipedia section on the composite Simpson's rule that covers the algorithm in detail.

I would also check out the SciPy implementation of Simpson's rule. Click here for the source code, which is well-documented.