MNUM - Simpson’s Rule in Python

April 11, 2008 in FEUP, Programming


def simpson(f, a, b, n):
    "Approximate the definite integral of f from a to b by Simpson's rule."

    if n % 2 != 0:
        print "Ups: n must be even!"
        return -1
        
    h  = (float(b) - a)/n
    
    si = 0.0
    sp = 0.0
    
    for i in range(1, n, 2):
        xk = a + i*h
        si += f(xk)
    
    for i in range(2, n, 2):
        xk = a + i*h
        sp += f(xk)
        
        
    s = 2*sp + 4*si + f(a) + f(b)

    return (h/3)*s

def f(x):
    return x**4

ni = 50
nf = 1000000
n = ni
a = -20
b = 0
s = []
qc = []
ec = []
t = 0.0
div = 0.0
i = 0.0

while n < nf:
    t = simpson(f, a, b, n)
    s.append(t)
    n = n * 2

for i in range(0, len(s)-2):
    div = (s[i+2] - s[i+1])
    if div == 0:
        break
    
    t = (s[i+1] - s[i]) / div
    qc.append(t)
    t = div / 15
    ec.append(t)

for i in range(0, len(qc)):
    print "%.12f, %.12f, %.12f => qc=%.12f, e=%.12f" % (s[i], s[i+1], s[i+2], qc[i], ec[i])

3 Replies

Tiago Boldt Pereira de Sousa
April 11, 2008 at 8:16 pm

Já vi que andamos no mesmo curso

Bem, já que todas as semanas fazes isso em python, hoje deixo a minha solução em C/C++. Só tens que mudar a libraria para imprimir os resultados. Precisa que seja definida a funcão f(x) e o seu integral If(x).

Resolve este exemplo instantaneamente em 2 ou 3 passos

#include

using namespace std;

long double f(long double x) { return x*x*x; }
long double If(long double x) { return x*x*x*x/4; }

long double simpson(double a, double b, int n){
long double h = (b-a)/n;
long double result=f(a)+f(b);

for(int steps=1; steps<n ; steps++)
result+=(steps%2)? 4*f(a+steps*h) : 2*f(a+steps*h);

return result*h/3;
}

int main(){
cout << “Integral de x^3: ” << If(20)-If(0) << endl;
cout << “Pela regra de Simpson: ” << simpson(0, 20, 10) << endl;
return 0;
}
lrei
April 11, 2008 at 9:30 pm

hehehe nice

I should probably justify the (apparent) complexity of my algorithm:

- The reason for the two cycles is that it removes the need for an “if” (or ?) making it faster.

- The reason for the multiplication (by two or four) to only be done in the end is to reduce errors

Both were instructions from the teacher, I take no credit for it.

As for the final part (which is outside of the simpson function) it’s all about errors and stuff. I could explain it but it’s beyond the scope of the post

1 MNUM - SimpsonCsv « LuisRei.com Pingback on May 10th, 2008 at 3:45 pm
[...] May 10, 2008 in FEUP, Programming This one instead reads the values from a CSV file containing experimental data. Link to my previous implementation of Simspson’s rule. [...]