Anonymous user
Numerical integration: Difference between revisions
GP
(GP) |
|||
Line 2,122:
let t = integrate square 0. 1. 10 trapezium
let s = integrate square 0. 1. 10 simpson</lang>
=={{header|PARI/GP}}==
Note also that double exponential integration is available as <code>intnum(x=a,b,f(x))</code> and Romberg integration is available as <code>intnumromb(x=a,b,f(x))</code>.
<lang parigp>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)</lang>
Results:
<pre>Testing function (x)->x^3 on [0, 1] with 100 intervals:
0.2450249999999999998
0.2499874999999999998
0.2550249999999999998
0.2500249999999999998
0.2499999999999999999
time = 0 ms.
Testing function (x)->1/x on [1, 100] with 1000 intervals:
4.654991057514676000
4.604762548678375026
4.556981057514676011
4.605986057514676146
4.605170384957142170
time = 15 ms.
Testing function (x)->x on [0, 5000] with 5000000 intervals:
12499997.49999919783
12499999.99999917123
12500002.49999919783
12499999.99999919783
12499999.99999923745
time = 29,141 ms.
Testing function (x)->x on [0, 6000] with 6000000 intervals:
17999996.99999869563
17999999.99999864542
18000002.99999869563
17999999.99999869563
17999999.99999863097
time = 34,820 ms.</pre>
=={{header|Pascal}}==
| |||