Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Double integration

Subject: Double integration

From: Marcio Barbalho

Date: 1 Dec, 2011 18:18:08

Message: 1 of 11

Dear friends

I've been trying to solve this double integration with the help of MatLab but apparently I am missing something.

INT (3,5) INT (0,SQRT(25-x^2)) SQRT((16x^2+16y^2-625)/(x^2+y^2-25)) dydx

The above is not in any standard syntax, but perhaps will communicate the actual problem. The integration on y is from 0 to SQRT(25-x^2) and the integration on x is from 3 to 5.

This is the m file:
%
syms x y
clc
firstint = int(sqrt((16*x*x+16*y*y-625)/(x*x+y*y-25)),y,0,sqrt(25-x^2))
answer = int(firstint,x,3,5)
%

and this is the error message I have got:

??? Attempt to call constructor double with incorrect letter case.
Error in ==> sym.sym>sym.sym/symnumeric at 152
                            S{k} = symr(double(x(k)));
Error in ==> sym.sym>sym.sym at 116
                S = cell2sym(S,symnumeric(x,a));
Error in ==> sym.times at 8
    A = sym(A);
Error in ==> sym.eval at 14
    s = evalin('caller',vectorize(map2mat(char(x))));

Many thanks

Subject: Double integration

From: Roger Stafford

Date: 1 Dec, 2011 22:26:08

Message: 2 of 11

"Marcio Barbalho" <marciobarbalho@yahoo.com> wrote in message <jb8gd0$52k$1@newscl01ah.mathworks.com>...
> firstint = int(sqrt((16*x*x+16*y*y-625)/(x*x+y*y-25)),y,0,sqrt(25-x^2))
- - - - - - - - - -
  I can't comment on your error messages, but I see that the form of your inner integral is that of an incomplete elliptic integral of the second kind, and this might be the cause of your difficulty performing symbolic integration.

  You can always try numerical integration techniques. However, even though the integrand is actually integrable, you will have to deal appropriately with the singularity that occurs along the periphery of the circle, x^2+y^2=5^2, presumably by a change of variables.

Roger Stafford

Subject: Double integration

From: Roger Stafford

Date: 2 Dec, 2011 03:05:10

Message: 3 of 11

"Marcio Barbalho" <marciobarbalho@yahoo.com> wrote in message <jb8gd0$52k$1@newscl01ah.mathworks.com>...
> INT (3,5) INT (0,SQRT(25-x^2)) SQRT((16x^2+16y^2-625)/(x^2+y^2-25)) dydx
- - - - - - - - - -
  In case you are interested in using numerical methods, in your problem you can change to polar coordinates which will yield an inner integral easily expressed in closed form, thereby avoiding any elliptic integrals. Using this expression you need only perform a single numerical integration with the outer integral having no singularities present.

  Changing to polar coordinates: x = r*cos(t) and y = r*sin(t) yields the double integral (using your notation):

INT (0,asec(5/3)) INT (3*sec(t),5) SQRT((625-16*r^2)/(25-r^2))*r dr dt

In the indefinite inner integral here, substituting u = sqrt(25-r^2) produces

 INT -SQRT(225+16*u^2) du

which has an easy closed form. This can be used in performing the outer integral numerically. (I don't see any closed form for this outer integration.)

Roger Stafford

Subject: Double integration

From: Marcio Barbalho

Date: 2 Dec, 2011 18:45:09

Message: 4 of 11

I've been trying with both Mupad and Maple... it just doesn't work. Of course I may be missing some tricks.
Apparent the result of the original double integration is 53.22. I would like to prove it.

Thank you

Subject: Double integration

From: Roger Stafford

Date: 2 Dec, 2011 19:55:09

Message: 5 of 11

"Marcio Barbalho" <marciobarbalho@live.com> wrote in message <jbb6bl$6af$1@newscl01ah.mathworks.com>...
> I've been trying with both Mupad and Maple... it just doesn't work. Of course I may be missing some tricks.
> Apparent the result of the original double integration is 53.22. I would like to prove it.
>
> Thank you
- - - - - - - - -
  It is my guess that you will never find a double integral solution to this particular problem of yours unless you do things numerically at some stage. I also doubt that Mupad or Maple will try experimenting with a polar coordinate change that could simplify the inner integral. What have you got against numerical methods, Marcio?

Roger Stafford

Subject: Double integration

From: Marcio Barbalho

Date: 2 Dec, 2011 20:42:09

Message: 6 of 11

"Roger Stafford" wrote
> - - - - - - - - -
> It is my guess that you will never find a double integral solution to this particular problem of yours unless you do things numerically at some stage. I also doubt that Mupad or Maple will try experimenting with a polar coordinate change that could simplify the inner integral. What have you got against numerical methods, Marcio?
>
> Roger Stafford

I have abso-bloody-lutely nothing against numerial integration. In actuality, I love numerical methods, of course, that does not mean I am an expert. I am just a junior engineer who loves maths and eventually find it interesting to play with calculus on rainy weekends.

Alright... my question is... how do I solve the inner integral numerically if the upper limit is not a number but an expression? I tried with the substitutions you proposed it pushed me to the drawing board, meaning it hit some point that neither mupad nor maple could go any further.

Thank you

Subject: Double integration

From: Roger Stafford

Date: 2 Dec, 2011 22:29:08

Message: 7 of 11

"Marcio Barbalho" <marciobarbalho@live.com> wrote in message <jbbd71$af$1@newscl01ah.mathworks.com>...
> I have abso-bloody-lutely nothing against numerial integration. In actuality, I love numerical methods, of course, that does not mean I am an expert. I am just a junior engineer who loves maths and eventually find it interesting to play with calculus on rainy weekends.
>
> Alright... my question is... how do I solve the inner integral numerically if the upper limit is not a number but an expression? I tried with the substitutions you proposed it pushed me to the drawing board, meaning it hit some point that neither mupad nor maple could go any further.
>
> Thank you
- - - - - - - - - -
  OK, do you agree with me that changing to polar coordinates gives you this equivalent double integral:

 int (0,asec(5/3)) int (3*sec(t),5) sqrt((625-16*r^2)/(25-r^2))*r dr dt ?

I'm not certain but I don't believe Mupad or Maple can do that for you. It takes a little calculus on your part.

  However, I strongly suspect that Maple can obtain a formula in terms of t for that inner integral taken with respect to r:

 int (3*sec(t),5) sqrt((625-16*r^2)/(25-r^2))*r dr

If I can do it using the substitution u = sqrt(25-r^2) as I showed you earlier, presumably Maple can also do it with some trick or other. If not, here is the basic form out of my integral table that will allow you to solve that indefinite integral in terms of u I mentioned:

 int sqrt(a^2+z^2) dz = 1/2*z*sqrt(a^2+z^2)+1/2*a^2*log(z+sqrt(a^2+z^2))

From that you can derive a solution to the above inner definite integral. Either way, this inner integral does not require numerical integration. The result is a specific formula in terms of t for the inner integral.

  Using the resulting formula for the inner definite integral as a function of t, you can then use one of matlab's numerical quadrature functions, 'quad', 'quadgk', or 'quadl', to do a single numerical (outer) integration with respect to t using the inner integral expression as an integrand between the two limits t=0 and t=asec(5/3). You do need to do a little study on these routines to be able to set up the desired error tolerances, etc.

  As I mentioned earlier there should be no trouble with singularities if you proceed along the above lines. On the other hand, if you were to integrate with respect to t first and then with respect to r afterwards, I believe you would have trouble with a singularity at r = 5, so it is advisable to integrate with respect to r first.

  Does all that make sense? If not, just ask for further help.

Roger Stafford

Subject: Double integration

From: Roger Stafford

Date: 4 Dec, 2011 01:27:09

Message: 8 of 11

"Marcio Barbalho" <marciobarbalho@live.com> wrote in message <jbb6bl$6af$1@newscl01ah.mathworks.com>...
> Apparent the result of the original double integration is 53.22. I would like to prove it.
- - - - - - - - -
  Marcio, I got around to applying those techniques I have described to you earlier in this thread and arrived at a somewhat higher value than the 53.22 you quoted. My answer was 53.573533745 . How did you arrive at your value?

  As I stated, there is a closed form for the inner integral using polar coordinates, and the outer integral was done numerically using one of matlab's quadrature routines with the inner integral formula as an integrand.

Roger Stafford

Subject: Double integration

From: Marcio Barbalho

Date: 4 Dec, 2011 02:03:09

Message: 9 of 11

"Roger Stafford" wrote in message <jbei9d$feq$1@newscl01ah.mathworks.com>...
> "Marcio Barbalho" <marciobarbalho@live.com> wrote in message <jbb6bl$6af$1@newscl01ah.mathworks.com>...
> > Apparent the result of the original double integration is 53.22. I would like to prove it.
> - - - - - - - - -
> Marcio, I got around to applying those techniques I have described to you earlier in this thread and arrived at a somewhat higher value than the 53.22 you quoted. My answer was 53.573533745 . How did you arrive at your value?
>
> As I stated, there is a closed form for the inner integral using polar coordinates, and the outer integral was done numerically using one of matlab's quadrature routines with the inner integral formula as an integrand.
>
> Roger Stafford

Very nice! I shall give it another go tomorrow. I got 53.22 using a calculator. No tricks, just typed the original double integral and got the result 2 min later. I set the calculator to work with 2 decimal places, otherwise it would take much longer to get to a 'better' result (like yours).

Thank you

Subject: Double integration

From: Steven_Lord

Date: 5 Dec, 2011 04:29:56

Message: 10 of 11



"Marcio Barbalho" <marciobarbalho@live.com> wrote in message
news:jbbd71$af$1@newscl01ah.mathworks.com...
> "Roger Stafford" wrote
>> - - - - - - - - -
>> It is my guess that you will never find a double integral solution to
>> this particular problem of yours unless you do things numerically at some
>> stage. I also doubt that Mupad or Maple will try experimenting with a
>> polar coordinate change that could simplify the inner integral. What
>> have you got against numerical methods, Marcio?
>>
>> Roger Stafford
>
> I have abso-bloody-lutely nothing against numerial integration. In
> actuality, I love numerical methods, of course, that does not mean I am an
> expert. I am just a junior engineer who loves maths and eventually find it
> interesting to play with calculus on rainy weekends.
>
> Alright... my question is... how do I solve the inner integral numerically
> if the upper limit is not a number but an expression?

http://www.mathworks.com/help/techdoc/ref/quad2d.html

--
Steve Lord
slord@mathworks.com
To contact Technical Support use the Contact Us link on
http://www.mathworks.com

Subject: Double integration

From: Roger Stafford

Date: 6 Dec, 2011 17:42:08

Message: 11 of 11

"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

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us