M242_lab31.html

M242 Monte Carlo 9/09

> with(Student[Calculus1]):

> f:=x-> 1/(1+x^2);

(1)

> T:= proc(N, a, b)
local h, k, s, exva;
h:= (b-a)/N:
s:= 2*sum(f(a+k*h), k=1..N-1):
exva:= int( f(x), x=a..b):
print( exva, evalf([exva, (f(a) + s + f(b))*h/2]));
end:

> T(100,0,1);

(2)

> ApproximateInt(f(x), x=0..1, method = trapezoid, partition=10); m:=evalf(%);

(3)

Practice Problems Today:

a) Try the case when f(x) = 1/x with the interval = [1,2].

b) Simpson's rule that we'll work on next lecture goes like this:

The interval is broken up into 2n (even) equal subsintervals. The weights distributed on the points

a= x[0], x[1], x[2], x[3], x[4], ... x[2n-1], x[2n] are respectively 1, 4, 2, 4, 2, ..., 4, 1.

In other words, the odd-indexed x-values are assigned weights 4, all the 'interior' even indexed x-values are assigned 2. The left and right end points are given weight = 1.

Simpson's Rule is: h:= (b-a)/(2*n), S:= (f(a) + 4*f(x1) + 2*f(x2) + 4*f(x3) + ... + 4*f(x_[n-1]) + f(b))*h/3.

Edit the procedure above, sum up the interior y-values by separating the indices into odd and even ones. Finally add them together with f(a) and f(b) as shown above.

Run the new proc on f(x) = 1/x to see if your answer agrees the one below.

> ApproximateInt(f(x), x=0..1, method = simpson, partition=10); m:=evalf(%);

(4)