Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Double integration
Date: Tue, 6 Dec 2011 17:42:08 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 22
Message-ID: <jblk5g$1lo$1@newscl01ah.mathworks.com>
References: <jb8gd0$52k$1@newscl01ah.mathworks.com> <jb9f96$ed2$1@newscl01ah.mathworks.com> <jbb6bl$6af$1@newscl01ah.mathworks.com> <jbbaet$jsd$1@newscl01ah.mathworks.com> <jbbd71$af$1@newscl01ah.mathworks.com> <jbhhc3$bau$1@newscl01ah.mathworks.com>
Reply-To: <HIDDEN>
NNTP-Posting-Host: www-00-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: newscl01ah.mathworks.com 1323193328 1720 172.30.248.45 (6 Dec 2011 17:42:08 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Tue, 6 Dec 2011 17:42:08 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1187260
Xref: news.mathworks.com comp.soft-sys.matlab:751573

"Steven_Lord" <slord@mathworks.com> wrote in message <jbhhc3$bau$1@newscl01ah.mathworks.com>...
> http://www.mathworks.com/help/techdoc/ref/quad2d.html
- - - - - - - -
  Hi Marcio and Steve.  In the case of that particular double integral I suggest that, rather than doing a fully two-dimensional numerical integration with 'quad2d', it may be preferable to transform to polar coordinates and take advantage of the fact that the resulting inner integrand now has a known analytic solution, so the problem can be reduced to a single numerical integration.

  With the change to polar coordinates, x = r*cos(t), y = r*sin(t), the problem assumes the form

int(int('k*sqrt(b^2-r^2)/sqrt(a^2-r^2)*r','r',h*sec(t),a),'t',0,asec(a/h))

where h = 3, k = 4, a = 5, and b = 6.25 .  Particularly if r is transformed to u = sqrt(a^2-r^2), a symbolic solution can readily be determined for this inner integrand.  When evaluated at the two limits of integration the resulting function, though no longer approaching infinity at the upper limit, nevertheless still possesses an infinite derivative as t approaches asec(a/h), something that may disturb a numerical integration procedure.  The transformation sin(v) = h/sqrt(a^2-h^2)*tan(t) remedies that also.

  The final result of all this is a single integral from v = 0 to v = pi/2 of the function

 f(v) = k*h*sqrt(a^2-h^2)/2*cos(v)./(a^2-(a^2-h^2)*cos(v).^2) .* ...
        (sqrt(a^2-h^2)*cos(v).*sqrt(b^2-a^2+(a^2-h^2)*cos(v).^2) ...
          + (b^2-a^2)*asinh(sqrt(a^2-h^2)/sqrt(b^2-a^2)*cos(v)));

This is a very well-behaved curve without singularities or infinite slopes.  In fact if plotted all the way from v = -pi/2 to +pi/2 it looks remarkably like a gaussian bell-shaped curve.  On my system the integral yields an answer of 53.57353374498558 when the tolerance is set to a tight value.

  The advantage this procedure has over a double numerical integration is presumably that far fewer points need be evaluated in the one-dimensional case, and there is no difficulty with the singularity that occurs all the way along the circular boundary.  On the other hand, there is admittedly more work for the user to do in the necessary algebra involved in performing these transformations.

Roger Stafford