Page 182 - Numerical methods for chemical engineering

168 4 Initial value problems

for k = 1:length(f)
var1 = x(k)ˆ2+yˆ2;
if(var1 > 1)
f(k) = 0;
else
f(k) = x(k)ˆ2 + 2*yˆ2 − 2*x(k)*y;
end
end
return;
The value of the integral is computed by

int 2D = dblquad(‘calc integrand 2D’,-1,1,-1,1),
int 2D = 2.3562

Monte Carlo integration
For integrations over very complex domains or in spaces of very high dimension, we use
Monte Carlo integration, a stochastic technique in which points are generated at random,
and contribute to the integral if they lie within the domain of integration. A more efficient
implementation of Monte Carlo integration, based on importance sampling, is discussed in
Chapter 7, and applications of this method to Bayesian statistical analysis are considered in
Chapter 8. Here, we present a simple Monte Carlo algorithm to compute (4.82). We define
the indicator function for the unit circle
2 2
1, x + y ≤ 1
r≤1 (x, y) = (4.83)
0, otherwise
2
2
Then, for f (x, y) = x + 2y − 2xy, the integral is
' 1 ' 1
I D = f (x, y) r≤1 (x, y)dxdy (4.84)
−1 −1
Using rand, a function that returns a uniformly-distributed random number in [0, 1], we
[k]
[k]
generate a long sequence of values (x , y ) and compute the average value over this
sequence of the integrand of (4.84),
1 [k] [k] [k] [k]
N s
f s = f x , y r≤1 x , y (4.85)
N s
k=1
For large N s , this sampled average should agree with the exact value
1 ' 1 ' 1
f exact = f (x, y) r≤1 (x, y)dxdy (4.86)
A sd −1 −1
where the area A sd of the sampled domain is

' 1 ' 1
A sd = (1)dxdy = 4 (4.87)
−1 −1