29 lines
767 B
Text
29 lines
767 B
Text
rectLeft(f, a, b, n)={
|
|
sum(i=0,n-1,f(a+(b-a)*i/n), 0.)*(b-a)/n
|
|
};
|
|
rectMid(f, a, b, n)={
|
|
sum(i=1,n,f(a+(b-a)*(i-.5)/n), 0.)*(b-a)/n
|
|
};
|
|
rectRight(f, a, b, n)={
|
|
sum(i=1,n,f(a+(b-a)*i/n), 0.)*(b-a)/n
|
|
};
|
|
trapezoidal(f, a, b, n)={
|
|
sum(i=1,n-1,f(a+(b-a)*i/n), f(a)/2+f(b)/2.)*(b-a)/n
|
|
};
|
|
Simpson(f, a, b, n)={
|
|
my(h=(b - a)/n, s);
|
|
s = 2*sum(i=1,n-1,
|
|
2*f(a + h * (i+1/2)) + f(a + h * i)
|
|
, 0.) + 4*f(a + h/2) + f(a) + f(b);
|
|
s * h / 6
|
|
};
|
|
test(f, a, b, n)={
|
|
my(v=[rectLeft, rectMid, rectRight, trapezoidal, Simpson]);
|
|
print("Testing function "f" on ",[a,b]," with "n" intervals:");
|
|
for(i=1,#v, print("\t"v[i](f, a, b, n)))
|
|
};
|
|
# \\ Turn on timer
|
|
test(x->x^3, 0, 1, 100)
|
|
test(x->1/x, 1, 100, 1000)
|
|
test(x->x, 0, 5000, 5000000)
|
|
test(x->x, 0, 6000, 6000000)
|