{"id":31,"date":"2006-04-26T09:56:43","date_gmt":"2006-04-26T14:56:43","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/?p=31"},"modified":"2018-01-08T16:13:17","modified_gmt":"2018-01-08T21:13:17","slug":"two-dimensional-integration-over-a-general-domain","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2006\/04\/26\/two-dimensional-integration-over-a-general-domain\/","title":{"rendered":"Two-Dimensional Integration over a General Domain"},"content":{"rendered":"<div xmlns:mwsh=\"https:\/\/www.mathworks.com\/namespace\/mcode\/v1\/syntaxhighlight.dtd\" class=\"content\">\r\n   <introduction>\r\n      <p>There have been several requests on the <a>MATLAB newsgroup<\/a> to calculate the area of a function on a non-rectangular domain.  MATLAB's function <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/ref\/dblquad.html\"><tt>dblquad<\/tt><\/a> uses the technique that I am going to illustrate today for rectangular domains. And I just noticed <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/loadFile.do?objectId=10639&amp;ref=rssfeed&amp;id=mostRecentFiles\">this entry<\/a> on the File Exchange for doing 2-dimensional quadrature over a non-rectangular domain.  It sounds similar but I have not\r\n         tried it out myself.\r\n      <\/p>\r\n   <\/introduction>\r\n   <h3>Contents<\/h3>\r\n   <div>\r\n      <ul>\r\n         <li><a href=\"#2\">Rectangular Area<\/a><\/li>\r\n         <li><a href=\"#3\">Rectangular Domain<\/a><\/li>\r\n         <li><a href=\"#5\">Elliptical Domain<\/a><\/li>\r\n         <li><a href=\"#11\">Summary<\/a><\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <p>Do a little housekeeping first.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">s = warning(<span style=\"color: #A020F0\">'off'<\/span>,<span style=\"color: #A020F0\">'MATLAB:quadl:MinStepSize'<\/span>);<\/pre><h3>Rectangular Area<a name=\"2\"><\/a><\/h3>\r\n   <p>To compute the area of a rectangle via integration, we can call <tt>dblquad<\/tt> as I mentioned earlier.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">a = 4;\r\nb = 3;\r\ndblquad(@(x,y) ones(size(x)), 0, a, 0, b)\r\n\r\n<span style=\"color: #228B22\">%<\/span>\r\n<span style=\"color: #228B22\">% and we see the result we expect: 3*4 = 12.<\/span><\/pre><pre style=\"font-style:oblique\">ans =\r\n    12\r\n<\/pre><h3>Rectangular Domain<a name=\"3\"><\/a><\/h3>\r\n   <p>To integrate something more interesting over that rectangular region, we simply supply a different function.  For example,\r\n      let's have <tt>x*y<\/tt> be the integrand.  We can do this analytically first.\r\n   <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/31\/quadN_eq394653.png\"> <\/p>\r\n   <p>Let's evaluate the integral numerically now with <tt>dblquad<\/tt> and compare the answer to the analytic one.  I am using the default tolerance and the <tt>quadl<\/tt> method for integration.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">tol = [];\r\nqnumeric = dblquad(@(x,y) x*y, 0, a, 0, b, tol, @quadl)\r\nqexact = (a^2 * b^2) \/ 4<\/pre><pre style=\"font-style:oblique\">qnumeric =\r\n    36\r\nqexact =\r\n    36\r\n<\/pre><h3>Elliptical Domain<a name=\"5\"><\/a><\/h3>\r\n   <p>The equation for an ellipse can be written as:<\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/31\/quadN_eq15613.png\"> <\/p>\r\n   <p>and we can rewrite <tt>y<\/tt> as a function of <tt>x<\/tt> for the upper-right quadrant as\r\n   <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/31\/quadN_eq14076.png\"> <\/p>\r\n   <p>We can now write the equation for the ellipse's area by integrating the identity integrand from 0 to the right limit for each\r\n      value of <tt>x<\/tt>.\r\n   <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/31\/quadN_eq96554.png\"> <\/p>\r\n   <p>and now we are ready to set up a function handle representing the inner integral.  Let's plot the function.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">fh = @(x) quadl(@(y) ones(size(y)), 0, b*sqrt(1-(x\/a).^2));\r\n\r\nx = 0:a\/100:a;\r\nell = arrayfun(fh, x);\r\nplot(x,ell); title(<span style=\"color: #A020F0\">'Ellipse boundary in first quadrant'<\/span>); axis <span style=\"color: #A020F0\">equal<\/span><\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/31\/quadN_01.png\"> \r\n<p>The reason I used <kbd>arrayfun<\/kbd> to evaluate <kbd>fh<\/kbd> above is the same reason why I can't just use <kbd>quadl<\/kbd> for the outer integral.  <kbd>fh<\/kbd> is not a vectorized function and I can't make it into one.  To evaluate <kbd>fh<\/kbd> for each value of <kbd>x<\/kbd>, I am calculating an integral at each point and the limit of integration is a function of <kbd>x<\/kbd>.  In addition, for each value of <kbd>x<\/kbd>, 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.\r\n<\/p>\r\n<p>Simply calling <tt>quadl<\/tt> with the function defined by <tt>fh<\/tt> 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 <tt>quadvec<\/tt>, the code I'm going to show you next, does (and it's similar to the inner workings of <tt>dblquad<\/tt> as well). \r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">dbtype <span style=\"color: #A020F0\">quadvec<\/span><\/pre><pre style=\"font-style:oblique\">\r\n1     function q = quadvec(f, varargin)\r\n2     % QUADVEC calls quad but assumes you can't vectorize your input function.\r\n3     %   QUADVEC has the same calling syntax as QUAD but doesn't require that\r\n4     %   the integrand be vectorized.  This makes QUADVEC suitable to perform\r\n5     %   double-, triple-, and higher order integrals over a non-rectangular\r\n6     %   domain.\r\n7     %   Q = QUADVEC(FUN,A,B) tries to approximate the integral of scalar-valued\r\n8     %   function FUN from A to B to within an error of 1.e-6 using recursive\r\n9     %   adaptive Simpson quadrature. FUN is a function handle. \r\n10    %\r\n11    %   Example:\r\n12    %   Compute area of ellipse by computing the result in one quadrant and\r\n13    %   then multiplying by 4.  Note that the inner integrand has integration\r\n14    %   limits that are a function of the semi-major and semi-minor ellipse\r\n15    %   axes, a and b.  \r\n16    %   First we set up the function handle to the inner integrand, fh.\r\n17    %   a = 4; b = 3;\r\n18    %   fh = @(x) quadl(@(y) ones(size(y)), 0, b*sqrt(1-(x\/a).^2));\r\n19    %   ea = 4 * quadvec(fh , 0, a);\r\n20    %   ea-pi*a*b  % compare integration result to closed-form solution\r\n21    \r\n22      q = quadl(@g, varargin{:}); % like quadl, but supplies g as the argument\r\n23      function y = g(X) % make f into a \"vectorized\" function\r\n24        y = zeros(size(X));\r\n25        for i = 1:numel(X)\r\n26          y(i) = f(X(i)); % this f refers to the argument of quadvec\r\n27        end\r\n28      end\r\n29    end\r\n\r\n<\/pre><p>You can see that we make sure we call the integrand for each value separately and we achieve this by making our function handle\r\n      \"vectorized\" by creating a function handle to a nested function, and using this new function, <tt>g<\/tt> as the integrand.  Looking at <tt>g<\/tt>, you can see, starting on line 25, the we use a for loop to iterate over all the input values.\r\n   <\/p>\r\n   <p>We're now ready to calculate the ellipse's area.  Remember, since we are only working in one quadrant, we need to multiply\r\n      the result by 4.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">ellipseArea = 4 * quadvec(fh,0,a)<\/pre><pre style=\"font-style:oblique\">ellipseArea =\r\n   37.6991\r\n<\/pre><p>and we compare this to the analytic formula:<\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/31\/quadN_eq2862.png\"> <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">ellipseAreaAnalytic = pi*a*b<\/pre><pre style=\"font-style:oblique\">ellipseAreaAnalytic =\r\n   37.6991\r\n<\/pre><p>Do a little final housekeeping.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">warning(s);<\/pre><h3>Summary<a name=\"11\"><\/a><\/h3>\r\n   <p>Today I showed you how to perform integration in higher dimensions in MATLAB with the introduction of <tt>quadvec<\/tt>.  You can see that <kbd>quadvec<\/kbd> ensures that we call the integrand one element at a time to allow us to cover\r\n      cases in which the integrand isn't vectorized. Note that you can use <tt>quadvec<\/tt> on an integrand in 1 dimension, but performance will probably be worse than using <tt>quadl<\/tt>, assuming the integrand is vectorized.\r\n   <\/p>\r\n   <p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br>\r\n      Published with MATLAB&reg; 7.2<br><\/p>\r\n<\/div>\r\n<!--\r\n##### SOURCE BEGIN #####\r\n%% Two-Dimensional Integration over a General Domain\r\n% There have been several requests on the <http:\/\/ MATLAB newsgroup>\r\n% to calculate the area of a function on a non-rectangular domain.  MATLAB's\r\n% function <https:\/\/www.mathworks.com\/help\/matlab\/ref\/dblquad.html |dblquad|>\r\n% uses the technique that I am going to illustrate today for rectangular\r\n% domains. And I just noticed\r\n% <https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/loadFile.do?objectId=10639&ref=rssfeed&id=mostRecentFiles this entry>\r\n% on the File Exchange for doing 2-dimensional quadrature over a\r\n% non-rectangular domain.  It sounds similar but I have not tried it out\r\n% myself.\r\n\r\n%%\r\n% Do a little housekeeping first.\r\ns = warning('off','MATLAB:quadl:MinStepSize');\r\n\r\n%% Rectangular Area\r\n% To compute the area of a rectangle via integration, we can call |dblquad|\r\n% as I mentioned earlier.\r\na = 4;\r\nb = 3;\r\ndblquad(@(x,y) ones(size(x)), 0, a, 0, b)\r\n\r\n%\r\n% and we see the result we expect: 3*4 = 12.\r\n\r\n%% Rectangular Domain\r\n% To integrate something more interesting over that rectangular region, we\r\n% simply supply a different function.  For example, let's have |x*y| be the\r\n% integrand.  We can do this analytically first.\r\n%\r\n% $$\\int_{0}^{b} \\int_{0}^{a} (x*y) dx dy = {x^{2}}|_{0}^{a} * {y^{2}}|_{0}^{b} \/ 4  = a^{2} * b^{2} \/ 4$$\r\n%\r\n\r\n%%\r\n% Let's evaluate the integral numerically now with |dblquad| and compare\r\n% the answer to the analytic one.  I am using the default tolerance and the\r\n% |quadl| method for integration.\r\ntol = [];\r\nqnumeric = dblquad(@(x,y) x*y, 0, a, 0, b, tol, @quadl)\r\nqexact = (a^2 * b^2) \/ 4\r\n\r\n%% Elliptical Domain\r\n% The equation for an ellipse can be written as:\r\n%\r\n% $$x^{2}\/a^{2} + y^{2}\/b^{2} = 1$$\r\n% \r\n% and we can rewrite |y| as a function of |x| for the upper-right quadrant\r\n% as\r\n%\r\n% $$y = b * \\surd(1 - (x\/a)^{2})$$\r\n% \r\n% We can now write the equation for the ellipse's area by integrating the\r\n% identity integrand from 0 to the right limit for each value of |x|.\r\n%\r\n% $$ \\int_{0}^{a} \\int_{0}^{b*\\surd (1-(x\/a)^{2})}  ( 1 ) dx dy $$\r\n%\r\n% and now we are ready to set up a function handle representing the inner\r\n% integral.  Let's plot the function.\r\n\r\nfh = @(x) quadl(@(y) ones(size(y)), 0, b*sqrt(1-(x\/a).^2));\r\n\r\nx = 0:a\/100:a;\r\nell = arrayfun(fh, x);\r\nplot(x,ell); title('Ellipse boundary in first quadrant'); axis equal\r\n\r\n%%\r\n% Simply calling |quadl| with the function defined by |fh| above won't work\r\n% because |quadl| expects to be able to evaluate the integrand in a\r\n% vectorized way and that can't be done with |fh|. Instead, we need an\r\n% integration method that doesn't assume the integrand is vectorized.  And\r\n% that's what |quadvec|, the code I'm going to show you next, does (and\r\n% it's similar to the inner workings of |dblquad| as well).   You can also\r\n% download |quadvec| <XXXXXX here>.\r\n \r\ndbtype quadvec\r\n\r\n%%\r\n% You can see that we make sure we call the integrand for each value\r\n% separately and we achieve this by making our function handle \"vectorized\"\r\n% by creating a function handle to a nested function, and using this new\r\n% function, |g| as the integrand.  Looking at |g|, you can see, starting on\r\n% line 25, the we use a for loop to iterate over all the input values.\r\n\r\n%%\r\n% We're now ready to calculate the ellipse's area.  Remember, since we are\r\n% only working in one quadrant, we need to multiply the result by 4.\r\nellipseArea = 4 * quadvec(fh,0,a)\r\n\r\n%%\r\n% and we compare this to the analytic formula:\r\n%\r\n% $$ \\pi * a * b $$\r\n%\r\nellipseAreaAnalytic = pi*a*b\r\n\r\n%%\r\n% Do a little final housekeeping.\r\nwarning(s);\r\n\r\n%% Summary\r\n% Today I showed you how to perform integration in higher dimensions in\r\n% MATLAB with the introduction of |quadvec|.  You can see that we used\r\n% a nested function to make sure we called the integrand one element at a\r\n% time to allow us to cover cases in which the integrand isn't vectorized.\r\n% Note that you can use |quadvec| on an integrand in 1\r\n% dimension, but performance will be decidedly worse than using |quadl|,\r\n% assuming the integrand is vectorized.\r\n\r\n##### SOURCE END #####\r\n-->","protected":false},"excerpt":{"rendered":"<p>\r\n   \r\n      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... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/loren\/2006\/04\/26\/two-dimensional-integration-over-a-general-domain\/\">read more >><\/a><\/p>","protected":false},"author":39,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[12],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/31"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/users\/39"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/comments?post=31"}],"version-history":[{"count":2,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/31\/revisions"}],"predecessor-version":[{"id":2610,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/31\/revisions\/2610"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/media?parent=31"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/categories?post=31"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/tags?post=31"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}