There have been several requests on the MATLAB newsgroup to calculate the area of a function on a non-rectangular domain. MATLAB's function dblquad uses the technique that I am going to illustrate today for rectangular domains. And I just noticed this entry on the File Exchange for doing 2-dimensional quadrature over a non-rectangular domain. It sounds similar but I have not tried it out myself.
Contents
Do a little housekeeping first.
s = warning('off','MATLAB:quadl:MinStepSize');
Rectangular Area
To compute the area of a rectangle via integration, we can call dblquad as I mentioned earlier.
a = 4; b = 3; dblquad(@(x,y) ones(size(x)), 0, a, 0, b) % % and we see the result we expect: 3*4 = 12.
ans =
12
Rectangular Domain
To integrate something more interesting over that rectangular region, we simply supply a different function. For example, let's have x*y be the integrand. We can do this analytically first.
Let's evaluate the integral numerically now with dblquad and compare the answer to the analytic one. I am using the default tolerance and the quadl method for integration.
tol = []; qnumeric = dblquad(@(x,y) x*y, 0, a, 0, b, tol, @quadl) qexact = (a^2 * b^2) / 4
qnumeric =
36
qexact =
36
Elliptical Domain
The equation for an ellipse can be written as:
and we can rewrite y as a function of x for the upper-right quadrant as
We can now write the equation for the ellipse's area by integrating the identity integrand from 0 to the right limit for each value of x.
and now we are ready to set up a function handle representing the inner integral. Let's plot the function.
fh = @(x) quadl(@(y) ones(size(y)), 0, b*sqrt(1-(x/a).^2)); x = 0:a/100:a; ell = arrayfun(fh, x); plot(x,ell); title('Ellipse boundary in first quadrant'); axis equal
The reason I used arrayfun to evaluate fh above is the same reason why I can't just use quadl for the outer integral. fh is not a vectorized function and I can't make it into one. To evaluate fh for each value of x, I am calculating an integral at each point and the limit of integration is a function of x. In addition, for each value of x, the integration method stops when the solution meets the correct numerical criteria and therefore each integration may not be performing the same number of function evaluations, nor will they necessarily evaluate the separate integrands at the same locations.
Simply calling quadl with the function defined by fh above won't work for the same reasons that I couldn't just vectorize my function evaluations of the integrand. Instead, we need an integration method that doesn't assume the integrand is vectorized. And that's what quadvec, the code I'm going to show you next, does (and it's similar to the inner workings of dblquad as well). You can also download quadvec here.
dbtype quadvec
1 function q = quadvec(f, varargin)
2 % QUADVEC calls quad but assumes you can't vectorize your input function.
3 % QUADVEC has the same calling syntax as QUAD but doesn't require that
4 % the integrand be vectorized. This makes QUADVEC suitable to perform
5 % double-, triple-, and higher order integrals over a non-rectangular
6 % domain.
7 % Q = QUADVEC(FUN,A,B) tries to approximate the integral of scalar-valued
8 % function FUN from A to B to within an error of 1.e-6 using recursive
9 % adaptive Simpson quadrature. FUN is a function handle.
10 %
11 % Example:
12 % Compute area of ellipse by computing the result in one quadrant and
13 % then multiplying by 4. Note that the inner integrand has integration
14 % limits that are a function of the semi-major and semi-minor ellipse
15 % axes, a and b.
16 % First we set up the function handle to the inner integrand, fh.
17 % a = 4; b = 3;
18 % fh = @(x) quadl(@(y) ones(size(y)), 0, b*sqrt(1-(x/a).^2));
19 % ea = 4 * quadvec(fh , 0, a);
20 % ea-pi*a*b % compare integration result to closed-form solution
21
22 q = quadl(@g, varargin{:}); % like quadl, but supplies g as the argument
23 function y = g(X) % make f into a "vectorized" function
24 y = zeros(size(X));
25 for i = 1:numel(X)
26 y(i) = f(X(i)); % this f refers to the argument of quadvec
27 end
28 end
29 end
You can see that we make sure we call the integrand for each value separately and we achieve this by making our function handle "vectorized" by creating a function handle to a nested function, and using this new function, g as the integrand. Looking at g, you can see, starting on line 25, the we use a for loop to iterate over all the input values.
We're now ready to calculate the ellipse's area. Remember, since we are only working in one quadrant, we need to multiply the result by 4.
ellipseArea = 4 * quadvec(fh,0,a)
ellipseArea = 37.6991
and we compare this to the analytic formula:
ellipseAreaAnalytic = pi*a*b
ellipseAreaAnalytic = 37.6991
Do a little final housekeeping.
warning(s);
Summary
Today I showed you how to perform integration in higher dimensions in MATLAB with the introduction of quadvec. You can see that quadvec ensures that we call the integrand one element at a time to allow us to cover cases in which the integrand isn't vectorized. Note that you can use quadvec on an integrand in 1 dimension, but performance will probably be worse than using quadl, assuming the integrand is vectorized.
Published with MATLAB® 7.2
good
Hi!
You are often making use of anonymous functions. Maybe you should state at some point that they are much slower than using a subfunction for the same task. Below I have put the code of a function that compares the cpu times in a simple example.
Markus
function testtime
L = 2000;
a = randn(L,10);
time0 = cputime;
for k=1:L
x = arrayfun(@(x)sin(x), a(k,:));
end
disp(sprintf(’Elapsed time is %g seconds.’, cputime-time0))
time0 = cputime;
for k=1:L
y = arrayfun(@subfun, a(k,:));
end
disp(sprintf(’Elapsed time is %g seconds.’, cputime-time0))
function y = subfun(x)
y = sin(x);
Performance of anonymous functions is sometimes less than for subfunctions currently. This is an area of performance (functions, nested functions, anonymous functions, and related) that we continue to work on. I find that for clariy and function signatures that show best the intent, nested and anonymous functions require less mental effort. The hope is that in future releases, the performance will continue to improve.
For some functions, we strive to eke out the maximum performance, but if we do so at the expense of significantly reduced clarity, we must think hard. There tradeoffs are subjective.
I’m a retired professor doing consulting on radar-raindrop relations. I have a choice of either resurrecting my decade-old Fortran programs or moving to a language like MatLab. My first hurdle in using Fortran is a call MSL to compute the surface area of a raindrop that is approximated as a triaxial ellipsoid. This allows me to explore possible normal modes of a raindrop having a given amount of oscillation energy that converts to surface energy at the oscillation end point of maximum area. Although I can find free Fortran compilers, the available IMSL seems to be unaffordable.
If I convert to MatLab, then I would use your quadvec in three dimensions. Any suggestions?
– Ken
triplequad ships as part of MATLAB. You might find that useful if you are integrating over fixed domains.
–Loren
Loren,
This function is really useful.
I have a little more complicated intergral problem:
f(x,y) = x+y;
I want to integrate f(x,y) with x varying from 0 to 1, and y varying from x to x+2;
Could you please give me some hint how to use the function?
This will solve a big problem for me!
Thank you very much.
Qi Chen-
Wouldn’t this do the trick?
fh = @(x) quadl(@(y) x+y, x, x+2); % plot y as a function of x x = 0:1/100:1; y = arrayfun(fh, x); plot(x,y); title('x vs. y'); triarea = quadvec(fh,0,1)–Loren
How do you multiply two anonymous functions with the same variable?The product has to go into the quad program.
Swadesh-
I am not sure what you mean but here’s a way to multiply to anonymous functions.
–Loren
Hi,
I tried to apply your function when the 2-d integral depends on a third variable, that is to say:
A=@(w) quadvec(fh,0,1);
In this case the “vectorized” function depends on w,
…
q =@(w)quadl(@(x)@g,varargin{:});
function y =g(x,w)
y = zeros(size(x));
for i = 1:numel(x)
y(i) = f{i}(x(i),w);(I wrote it like a cell array)
end
…
Unfortunately it doesn’t work…I can’t define q in an efficient way. May you give me any suggestion please?
Dear Loren,
I am not sure if this is a right place to post my question; however, my question is somehow related to this post.
In my study, I need to numerically evaluate a quadruple integral. However, Matlab only offers built in functions up to triple integral. Can you offer me some suggestion, please?
David-
Simply use the technique up above calling quadvec multiple times (i.e., have quadvec call quadvec, etc. up to your actual integrand).
–Loren
Loren,
How can I numerically integrate the function:
exp(1/(|x|^2 - 1), where |x| is the length of x, and x is within R^N for N=1 and N=2
I have been scratching my head trying to solve this problem and any help would be appreciated. Thanks
hi! Loren,
In this post, theres a small typo:
fh = @(x) quadl(@(y) x+y, x, x+2);
% plot y as a function of x
x = 0:1/100:1;
y = arrayfun(fh, x);
plot(x,y); title(’x vs. y’);
triarea = quadvec(fh,0,1)
I think the last line should be quadv(fh,0,1)
Please let me know.
Thanks
quadv does work. It was released after the original post. The quadvec code is available from the post itself. That last call is valid if you have made a copy of quadvec.
–Loren
Sorry Loren. I should have paid more attention. Thank you.
Hello,
I am M.Res student in Yangon University.My thesis title is 2 dimension elliptic equations.I have one problem about the matlab program.My equation has constant no(a(x))and (u’(x,y))’.I cannot get a correct 2 dimension graph.I have boundary condition for constant “a”.N=4.h=1/N+1((0,1)*(0,1)),for j=1:(N+1)/2…..alpha(j)=1;
for j=(N+1)/2+1:N+1…..alpha(j)=0.5;my equation is a’u'(x,y)+a(u’(x,y))’…..
I hope you can give ur suggestion for me.
I am waiting for your reply.
your sincerely
Thiri
Hi Loren,
Is there a way to perform double integration faster than using int(int(…)), where the limits of each integrand are different. I’m trying to perform double integration for a large array (2000×70), so by using the int(int(…)) approach, it is very time-consuming (~25min).
Hope you can help.
Best regards,
Kelly
Kelly-
int is for symbolic functions. You seem to have numeric values. I think there’s enough detail missing from your question that I can’t honestly help. I suggest you contact technical support with the entire M-file details.
–Loren
Loren,
Is there a way to carry out integration without using the function commands?
Thanks
Shawn
Shawn-
I don’t understand your question. I don’t know what a function command is.
–Loren
Dear Loren,
I want to numerically integrate a profit function with several integrales with the dblquad-command, therefore I ran two loops, and put the single variables from the loop as indices. The region I want to integrate is non-rectangular, that´s how I came to the quadvec-command. In the end I want to receive a big matrix of possible profits, where each element symbolizes a possible profit-outcome with the variables used from the loops. The variables from the loops are all integer and positives, but Matlab gives me the message: Subscript indices must either be real positive integers or logicals. I don´t really understand why. Therefore I tried to run my program with the quadvec-command which you posted here. But it gave the same error message. Do you have an idea, where I could have made a mistake?
Thanks in advance
Angela
Angela-
I have no idea. You should use a debugger with your code to see if that helps. And if you are still stuck, you might work with technical support to sort out the issue.
–Loren
Dear Loren,
I found some problems while using the quad function. For example, if I want to evaluate the expectation value of an exponential distributed variable, I approximately use the following:
>> quadl(@(x) x.*exp(-x), 0, 100)
ans =
1.0000
It works fine. However, 100 is just an approximation to infinity. So I wonder if I want a more precise answer, I should increase the upper limit. So I decide to use (0,1000). The answer is quite unexpected:
>> quadl(@(x) x.*exp(-x), 0, 1000)
ans =
1.9164e-036
Why does this happen?
The same thing also happens when I use dblquad function:
>> dblquad(@(x, y) y.*exp(-x).*exp(-y).*(y>a*x), 0, 100, 0, 100)
ans =
0.7500
>> dblquad(@(x, y) y.*exp(-x).*exp(-y).*(y>a*x), 0, 1000, 0, 1000)
ans =
4.5183e-025
Vincent-
I suspect it’s because the integration intervals widen when the integration limits widen, if the function is fairly flat, but might have some “peaks”, the peaks can be undersampled.
You might to contact technical support for more detail.
–Loren