{"id":825,"date":"2013-12-26T13:29:25","date_gmt":"2013-12-26T18:29:25","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/?p=825"},"modified":"2013-12-26T13:29:25","modified_gmt":"2013-12-26T18:29:25","slug":"double-integration-in-matlab-understanding-tolerances","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2013\/12\/26\/double-integration-in-matlab-understanding-tolerances\/","title":{"rendered":"Double Integration in MATLAB &#8211; Understanding Tolerances"},"content":{"rendered":"\r\n<div class=\"content\"><!--introduction--><p>In today's post, I am joined by Mike Hosea, a developer who occasionally works on integration routines for MATLAB.  In recent releases, we added new integration routines, including <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/ref\/integral2.html\"><tt>integral2<\/tt><\/a> for double integrals.  Today we'll talk about some nuances in using this routine and adjusting the absolute and relative tolerances.<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#4a1db13b-36aa-4d3c-9ca8-6fc2733369be\">Requirement for Calling Integration Functions<\/a><\/li><li><a href=\"#b6b25816-da83-4f37-97e4-24422f72ea3c\">Exploring Tolerances<\/a><\/li><li><a href=\"#c68026db-67da-443a-8ef9-763b7d05a9e3\">Adjust the Tolerance(s)<\/a><\/li><li><a href=\"#1442d89b-8f82-44a8-aab3-808a679ef87c\">Have You Been Successful with the Newer Integration Routines?<\/a><\/li><\/ul><\/div><h4>Requirement for Calling Integration Functions<a name=\"4a1db13b-36aa-4d3c-9ca8-6fc2733369be\"><\/a><\/h4><p>There is at least one requirement on the integrand when you call any integration function, and that is that it can be called with multiple inputs and it will return corresponding outputs, 1-for-1, with the inputs. So, to integrate $x * y$, you need to supply the <b>vectorized<\/b> form, i.e., <tt>x .* y<\/tt>.<\/p><p>Instead of creating a file with the integrand definition, we instead create an <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/matlab_prog\/anonymous-functions.html\">anonymous function<\/a>.  If you are unfamiliar with anonymous functions and want to learn more about them, please check the categories to the right of the blog for discussions and examples.<\/p><pre class=\"codeinput\">F0 = @(x,y) x .* y;\r\n<\/pre><p>We can now integrate <tt>F<\/tt> between 0 and 1 for both directions, <tt>x<\/tt> and <tt>y<\/tt>.<\/p><pre class=\"codeinput\">Q = integral2(F0,0,1,0,1)\r\n<\/pre><pre class=\"codeoutput\">Q =\r\n   0.250000000000000\r\n<\/pre><h4>Exploring Tolerances<a name=\"b6b25816-da83-4f37-97e4-24422f72ea3c\"><\/a><\/h4><p>Here is another example. First let's integrate the function <tt>F1<\/tt><\/p><pre class=\"codeinput\">F = @(x,y)exp(10*x).*tan(y)\r\n<\/pre><pre class=\"codeoutput\">F = \r\n    @(x,y)exp(10*x).*tan(y)\r\n<\/pre><p>between 0 and 1 for <tt>x<\/tt> and 0 and $\\pi\/4$ for <tt>y<\/tt>.<\/p><pre class=\"codeinput\">Q1 = integral2(F,0,1,0,pi\/4)\r\n<\/pre><pre class=\"codeoutput\">Q1 =\r\n     7.633444730985692e+02\r\n<\/pre><p>The analytic answer for this integral is<\/p><pre class=\"codeinput\">Q0 = log(sqrt(2))\/10*(exp(10) - 1)\r\n<\/pre><pre class=\"codeoutput\">Q0 =\r\n     7.633444758094897e+02\r\n<\/pre><p>and you can see that <tt>Q1<\/tt> differs from <tt>Q0<\/tt>, perhaps more than we might want.  We first adjust the numeric display to see enough decimal digits, before displaying the abolute error of the integration <tt>Q1<\/tt>.<\/p><pre class=\"codeinput\">format <span class=\"string\">long<\/span>\r\nintabserr1 = abs(Q1-Q0)\r\n<\/pre><pre class=\"codeoutput\">intabserr1 =\r\n     2.710920512072335e-06\r\n<\/pre><h4>Adjust the Tolerance(s)<a name=\"c68026db-67da-443a-8ef9-763b7d05a9e3\"><\/a><\/h4><p>One of the advantages of the <tt>integral2<\/tt> routine (and its companions for other dimensions) over most of the older routines is that they provide mixed relative and absolute error control, which is to say, you can supply two tolerances instead of just one. Not only can you, you should. You actually need both, or rather sometimes the one and sometimes the other, but it is good practice to think about what you want for both of them and to be aware of the default values in case you choose not to set them. What you would rarely want to do is change only the absolute error tolerance, <tt>AbsTol<\/tt>.<\/p><p>Now suppose we want more accuracy and try to crank down the absolute tolerance.<\/p><pre class=\"codeinput\">Q2 = integral2(F,0,1,0,pi\/4,<span class=\"string\">'AbsTol'<\/span>,1e-14)\r\nintabserr2 = abs(Q2 - Q0)\r\n<\/pre><pre class=\"codeoutput\">Q2 =\r\n     7.633444730985692e+02\r\nintabserr2 =\r\n     2.710920512072335e-06\r\n<\/pre><p>We got the same answer as before! The reason is that <tt>integral2<\/tt> uses a mixture of relative and absolute error control, and it is always the least restrictive of the two that matters. The default <tt>AbsTol<\/tt> is 1e-10, and the default <tt>RelTol<\/tt> is 1e-6. The absolute error is bigger than the absolute tolerance, true, but the relative error is<\/p><pre class=\"codeinput\">intrelerr2 = abs(Q2 - Q0)\/abs(Q0)\r\n<\/pre><pre class=\"codeoutput\">intrelerr2 =\r\n     3.551372411777180e-09\r\n<\/pre><p>and that's smaller than <tt>RelTol<\/tt>, so <tt>integral2<\/tt> returned that answer. Since the relative error is between 1e-8 and 1e-9, there should be about 8 correct significant digits, and indeed, we have that.<\/p><pre class=\"codeinput\">disp([Q0;Q2])\r\n<\/pre><pre class=\"codeoutput\">   1.0e+02 *\r\n   7.633444758094897\r\n   7.633444730985692\r\n<\/pre><p>The default <tt>RelTol<\/tt> is 1e-6, so let's try making that smaller, instead of <tt>AbsTol<\/tt>, to get more accuracy.<\/p><pre class=\"codeinput\">Q3 = integral2(F,0,1,0,pi\/4,<span class=\"string\">'RelTol'<\/span>,1e-12)\r\nintrelerr3 = abs(Q3 - Q0)\/abs(Q0)\r\n<\/pre><pre class=\"codeoutput\">Q3 =\r\n     7.633444758094937e+02\r\nintrelerr3 =\r\n     5.212639177138190e-15\r\n<\/pre><p>Now since the relative error is between 1e-14 and 1e-15, we should have about 14 significant digits correct, and we do<\/p><pre class=\"codeinput\">disp([Q0;Q3])\r\n<\/pre><pre class=\"codeoutput\">   1.0e+02 *\r\n   7.633444758094897\r\n   7.633444758094937\r\n<\/pre><p>Usually <tt>RelTol<\/tt> is the more important tolerance to manipulate, and the more useful one, because it indirectly controls the number of correct significant digits. But there is always a catch. Suppose the correct answer is zero. Then relative error is undefined! That's why <tt>integral2<\/tt> has to use the least restrictive of <tt>AbsTol<\/tt> and <tt>RelTol<\/tt>. The actual error condition it tries to meet is<\/p><pre>     ErrEst &lt;= max(AbsTol,RelTol*abs(Q))<\/pre><p>where <tt>ErrEst<\/tt> is an approximate upper bound on the absolute error.<\/p><p>It turns out that the ratio <tt>R = AbsTol\/RelTol<\/tt> is a cutoff. When <tt>abs(Q) &lt; R<\/tt>, <tt>AbsTol<\/tt> is the tolerance that matters. When <tt>abs(Q) &gt; R<\/tt>, <tt>RelTol<\/tt> is the one that matters. So it is good practice to set both tolerances to something reasonable by asking yourself two questions:<\/p><div><ol><li>If the true solution is small in magnitude, how many decimal places do I care about? The answer to that question tells you what <tt>AbsTol<\/tt> should be. If the answer is 10 places, then make <tt>AbsTol<\/tt> 1e-10 or perhaps 5e-11, so that the last digit should be correct when rounding to that point.<\/li><li>Assuming the answer is not too small, how many correct significant digits do I need? The answer to that question tells you what <tt>RelTol<\/tt> should be. Here again, if you would like 10 significant digits, make <tt>RelTol<\/tt> 1e-10 or perhaps 5e-11.<\/li><\/ol><\/div><p>So when you choose your own tolerances, the calls should look something like<\/p><pre class=\"codeinput\">Q4 = integral2(F,0,1,0,pi\/4,<span class=\"string\">'AbsTol'<\/span>,1e-12,<span class=\"string\">'RelTol'<\/span>,1e-10)\r\n<\/pre><pre class=\"codeoutput\">Q4 =\r\n     7.633444758095101e+02\r\n<\/pre><h4>Have You Been Successful with the Newer Integration Routines?<a name=\"1442d89b-8f82-44a8-aab3-808a679ef87c\"><\/a><\/h4><p>In this post, we've only covered the topic of setting tolerances to insure a reasonable estimate for your integrals.  We plan to talk more about which method to choose, order of integration (e.g., <tt>x<\/tt> then <tt>y<\/tt>), and dealing with singularities and discontinuities.<\/p><p>In the meantime, do you have any experiences to share about performing your own numerical integrations with MATLAB?  Let us know <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=825#respond\">here<\/a>.<\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_3f57334536714dd1ac4f9cbe71718f73() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='3f57334536714dd1ac4f9cbe71718f73 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 3f57334536714dd1ac4f9cbe71718f73';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        copyright = 'Copyright 2013 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }   \r\n     --> <\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_3f57334536714dd1ac4f9cbe71718f73()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n      the MATLAB code <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; R2013b<br><\/p><p class=\"footer\"><br>\r\n      Published with MATLAB&reg; R2013b<br><\/p><\/div><!--\r\n3f57334536714dd1ac4f9cbe71718f73 ##### SOURCE BEGIN #####\r\n%% Double Integration in MATLAB - Understanding Tolerances\r\n% In today's post, I am joined by Mike Hosea, a developer who occasionally\r\n% works on integration routines for MATLAB.  In recent releases, we added\r\n% new integration routines, including\r\n% <https:\/\/www.mathworks.com\/help\/matlab\/ref\/integral2.html |integral2|> for\r\n% double integrals.  Today we'll talk about some nuances in using this\r\n% routine and adjusting the absolute and relative tolerances.\r\n\r\n%% Requirement for Calling Integration Functions\r\n% There is at least one requirement on the integrand when you call any\r\n% integration function, and that is that it can be called with multiple\r\n% inputs and it will return corresponding outputs, 1-for-1, with the\r\n% inputs. So, to integrate $x * y$, you need to supply the *vectorized*\r\n% form, i.e., |x .* y|.\r\n\r\n%%\r\n% Instead of creating a file with the integrand definition, we instead\r\n% create an\r\n% <https:\/\/www.mathworks.com\/help\/matlab\/matlab_prog\/anonymous-functions.html\r\n% anonymous function>.  If you are unfamiliar with anonymous functions and\r\n% want to learn more about them, please check the categories to the right\r\n% of the blog for discussions and examples.\r\nF0 = @(x,y) x .* y;\r\n\r\n%%\r\n% We can now integrate |F| between 0 and 1 for both directions, |x| and\r\n% |y|.\r\nQ = integral2(F0,0,1,0,1)\r\n%% Exploring Tolerances\r\n% Here is another example. First let's integrate the function |F1|\r\nF = @(x,y)exp(10*x).*tan(y)\r\n\r\n%%\r\n% between 0 and 1 for |x| and 0 and $\\pi\/4$ for |y|.\r\nQ1 = integral2(F,0,1,0,pi\/4)\r\n\r\n%%\r\n% The analytic answer for this integral is\r\nQ0 = log(sqrt(2))\/10*(exp(10) - 1)\r\n\r\n%% \r\n% and you can see that |Q1| differs from |Q0|, perhaps more than we might\r\n% want.  We first adjust the numeric display to see enough decimal digits,\r\n% before displaying the abolute error of the integration |Q1|.\r\nformat long\r\nintabserr1 = abs(Q1-Q0)\r\n\r\n%% Adjust the Tolerance(s)\r\n% One of the advantages of the |integral2| routine (and its companions for\r\n% other dimensions) over most of the older routines is that they provide\r\n% mixed relative and absolute error control, which is to say, you can\r\n% supply two tolerances instead of just one. Not only can you, you should.\r\n% You actually need both, or rather sometimes the one and sometimes the\r\n% other, but it is good practice to think about what you want for both of\r\n% them and to be aware of the default values in case you choose not to set\r\n% them. What you would rarely want to do is change only the absolute error\r\n% tolerance, |AbsTol|.  \r\n%\r\n% Now suppose we want more accuracy and try to crank down the absolute\r\n% tolerance.\r\nQ2 = integral2(F,0,1,0,pi\/4,'AbsTol',1e-14)\r\nintabserr2 = abs(Q2 - Q0)\r\n\r\n%%\r\n% We got the same answer as before! The reason is that |integral2| uses a\r\n% mixture of relative and absolute error control, and it is always the\r\n% least restrictive of the two that matters. The default |AbsTol| is 1e-10,\r\n% and the default |RelTol| is 1e-6. The absolute error is bigger than the\r\n% absolute tolerance, true, but the relative error is\r\nintrelerr2 = abs(Q2 - Q0)\/abs(Q0)\r\n\r\n%%\r\n% and that's smaller than |RelTol|, so |integral2| returned that answer.\r\n% Since the relative error is between 1e-8 and 1e-9, there should be about\r\n% 8 correct significant digits, and indeed, we have that.\r\ndisp([Q0;Q2])\r\n\r\n%%\r\n% The default |RelTol| is 1e-6, so let's try making that smaller, instead\r\n% of |AbsTol|, to get more accuracy.\r\nQ3 = integral2(F,0,1,0,pi\/4,'RelTol',1e-12)\r\nintrelerr3 = abs(Q3 - Q0)\/abs(Q0)\r\n\r\n%%\r\n% Now since the relative error is between 1e-14 and 1e-15, we should have\r\n% about 14 significant digits correct, and we do\r\ndisp([Q0;Q3])\r\n\r\n%%\r\n% Usually |RelTol| is the more important tolerance to manipulate, and the\r\n% more useful one, because it indirectly controls the number of correct\r\n% significant digits. But there is always a catch. Suppose the correct\r\n% answer is zero. Then relative error is undefined! That's why |integral2|\r\n% has to use the least restrictive of |AbsTol| and |RelTol|. The actual\r\n% error condition it tries to meet is\r\n%\r\n%       ErrEst <= max(AbsTol,RelTol*abs(Q))\r\n%\r\n% where |ErrEst| is an approximate upper bound on the absolute error.\r\n%%\r\n% It turns out that the ratio |R = AbsTol\/RelTol| is a cutoff. When\r\n% |abs(Q) < R|, |AbsTol| is the tolerance that matters. When |abs(Q) > R|,\r\n% |RelTol| is the one that matters. So it is good practice to set both\r\n% tolerances to something reasonable by asking yourself two questions:\r\n%\r\n% # If the true solution is small in magnitude, how many decimal places do\r\n% I care about? The answer to that question tells you what |AbsTol| should\r\n% be. If the answer is 10 places, then make |AbsTol| 1e-10 or perhaps 5e-11,\r\n% so that the last digit should be correct when rounding to that point.\r\n% # Assuming the answer is not too small, how many correct significant\r\n% digits do I need? The answer to that question tells you what |RelTol|\r\n% should be. Here again, if you would like 10 significant digits, make\r\n% |RelTol| 1e-10 or perhaps 5e-11.\r\n%\r\n% So when you choose your own tolerances, the calls should look something\r\n% like\r\n\r\nQ4 = integral2(F,0,1,0,pi\/4,'AbsTol',1e-12,'RelTol',1e-10)\r\n\r\n%% Have You Been Successful with the Newer Integration Routines?\r\n% In this post, we've only covered the topic of setting tolerances to\r\n% insure a reasonable estimate for your integrals.  We plan to talk more\r\n% about which method to choose, order of integration (e.g., |x| then |y|),\r\n% and dealing with singularities and discontinuities.\r\n%\r\n% In the meantime, do you have any experiences to share about performing\r\n% your own numerical integrations with MATLAB?  Let us know\r\n% <https:\/\/blogs.mathworks.com\/loren\/?p=825#respond here>.\r\n\r\n\r\n##### SOURCE END ##### 3f57334536714dd1ac4f9cbe71718f73\r\n-->","protected":false},"excerpt":{"rendered":"<!--introduction--><p>In today's post, I am joined by Mike Hosea, a developer who occasionally works on integration routines for MATLAB.  In recent releases, we added new integration routines, including <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/ref\/integral2.html\"><tt>integral2<\/tt><\/a> for double integrals.  Today we'll talk about some nuances in using this routine and adjusting the absolute and relative tolerances.... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/loren\/2013\/12\/26\/double-integration-in-matlab-understanding-tolerances\/\">read more >><\/a><\/p>","protected":false},"author":39,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[59,26],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/825"}],"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=825"}],"version-history":[{"count":3,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/825\/revisions"}],"predecessor-version":[{"id":828,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/825\/revisions\/828"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/media?parent=825"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/categories?post=825"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/tags?post=825"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}