Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Double Integral Singularities quad2d
Date: Wed, 26 Oct 2011 15:01:14 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 235
Message-ID: <j897bq$e3s$1@newscl01ah.mathworks.com>
References: <j838vr$9nu$1@newscl01ah.mathworks.com> <41a15ffe-e999-41b9-8261-bd74f589cb23@l12g2000vby.googlegroups.com> <77e489a5-70a9-42a2-a577-3db5271345bc@eh5g2000vbb.googlegroups.com>
Reply-To: <HIDDEN>
NNTP-Posting-Host: www-01-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: newscl01ah.mathworks.com 1319641274 14460 172.30.248.46 (26 Oct 2011 15:01:14 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Wed, 26 Oct 2011 15:01:14 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 871220
Xref: news.mathworks.com comp.soft-sys.matlab:747254

Torsten <Torsten.Hennig@umsicht.fraunhofer.de> wrote in message <77e489a5-70a9-42a2-a577-3db5271345bc@eh5g2000vbb.googlegroups.com>...
> On 26 Okt., 09:11, "S " <simaher2...@yahoo.co.uk> wrote:
> > Torsten <Torsten.Hen...@umsicht.fraunhofer.de> wrote in message <05f56d39-0191-45b8-9fd8-286f5a444...@r28g2000yqj.googlegroups.com>...
> > > On 25 Okt., 14:13, "S " <simaher2...@yahoo.co.uk> wrote:
> > > > Torsten <Torsten.Hen...@umsicht.fraunhofer.de> wrote in message <d15f286c-7e52-4d35-8bad-9d77923e7...@p16g2000yqj.googlegroups.com>...
> > > > > On 25 Okt., 11:06, "S " <simaher2...@yahoo.co.uk> wrote:
> > > > > > Torsten <Torsten.Hen...@umsicht.fraunhofer.de> wrote in message <272cde32-3068-4756-baba-6252dcd80...@c1g2000vbw.googlegroups.com>...
> > > > > > > On 24 Okt., 17:18, "S " <simaher2...@yahoo.co.uk> wrote:
> > > > > > > > Torsten <Torsten.Hen...@umsicht.fraunhofer.de> wrote in message <36eea4f4-4887-4557-8c14-744e0c5bb...@q13g2000vbd.googlegroups.com>...
> > > > > > > > > On 24 Okt., 11:51, "S " <simaher2...@yahoo.co.uk> wrote:
> > > > > > > > > > Torsten <Torsten.Hen...@umsicht.fraunhofer.de> wrote in message <db464627-33b0-456e-907c-47f521590...@f36g2000vbm.googlegroups.com>...
> > > > > > > > > > > On 24 Okt., 11:24, "S " <simaher2...@yahoo.co.uk> wrote:
> > > > > > > > > > > > Torsten <Torsten.Hen...@umsicht.fraunhofer.de> wrote in message <41a15ffe-e999-41b9-8261-bd74f589c...@l12g2000vby.googlegroups.com>...
> > > > > > > > > > > > > On 24 Okt., 10:52, "S " <simaher2...@yahoo.co.uk> wrote:
> > > > > > > > > > > > > > Hi,
> >
> > > > > > > > > > > > > > I have a problem in solving the double integral below in matlab.
> >
> > > > > > > > > > > > > > As I increase f and g to larger values then matlab complains of singularity and unsuccessful integration.
> >
> > > > > > > > > > > > > > 3 - 2*cos(y).*cos(f*x+g*y) - cos((f-1)*x + g*y)  ./  4 - 2*cos(y).*(cos(y) + cos(x)) dxdy
> >
> > > > > > > > > > > > > > Over limits -pi to pi for x and -pi to pi for y.
> >
> > > > > > > > > > > > > > I have been using:
> > > > > > > > > > > > > > quad2d(@(x,y)my_func(x,y,f,g),-pi,pi,-pi,pi)
> >
> > > > > > > > > > > > > > Can anyone please help? How can I overcome this problem as I would like to integrate this function for various 'f' and 'g' up to ~2000.
> >
> > > > > > > > > > > > > > Thanks!
> >
> > > > > > > > > > > > > Since the period of your function becomes smaller and smaller with
> > > > > > > > > > > > > increasing values for f and g,
> > > > > > > > > > > > > you should integrate it analytically with the help of MATLAB's int.
> >
> > > > > > > > > > > > > Best wishes
> > > > > > > > > > > > > Torsten.
> >
> > > > > > > > > > > > Hi Torsten,
> >
> > > > > > > > > > > > Thanks for your reply.
> >
> > > > > > > > > > > > I take your point but I would like to integrate this numerically. Ideally, I would like to compute the integral for various f and g.
> >
> > > > > > > > > > > > Any ideas how I can do this?- Zitierten Text ausblenden -
> >
> > > > > > > > > > > > - Zitierten Text anzeigen -
> >
> > > > > > > > > > > Can you prescribe an initial x- and y-grid for MATLAB's quad2d ?
> > > > > > > > > > > If so, choose the increment in x and y very small to ensure that your
> > > > > > > > > > > highly oscillating function
> > > > > > > > > > > is sufficiently resolved.
> > > > > > > > > > > To get an impression of your function, you should plot it for high
> > > > > > > > > > > values of f and g.
> > > > > > > > > > > Do you see now why quad2d has difficulties to integrate it
> > > > > > > > > > > numerically ?
> >
> > > > > > > > > > > I repeat: Integrate your function _analytically_ using MATLAB's 'int'.
> >
> > > > > > > > > > > Best wishes
> > > > > > > > > > > Torsten.
> >
> > > > > > > > > > Thanks for your reply and so fast.
> >
> > > > > > > > > > However, I do not have the symbolic toolbox installed :-(
> >
> > > > > > > > > > Surely there must be another way?- Zitierten Text ausblenden -
> >
> > > > > > > > > > - Zitierten Text anzeigen -
> >
> > > > > > > > > Using
> > > > > > > > >http://integrals.wolfram.com
> > > > > > > > > ,I get
> > > > > > > > > 8*Pi^2 + sin(Pi*f)*sin(Pi*g)*(1/(g*(f-1)) + 8*g/(f*(1-g^2)))
> > > > > > > > > for your integral.
> > > > > > > > > But you should double-check the result.
> >
> > > > > > > > > Best wishes
> > > > > > > > > Torsten.
> >
> > > > > > > > Thanks for reply.
> >
> > > > > > > > @ Torsten. Very good website btw. However, the website is unable to integrate the function.- Zitierten Text ausblenden -
> >
> > > > > > > > - Zitierten Text anzeigen -
> >
> > > > > > > Use the website as follows:
> > > > > > > First integrate your function with respect to x by treating y as
> > > > > > > constant:
> > > > > > > int(3) = 3*x
> > > > > > > int(- 2*cos(y)*cos(f*x+g*y))=-2*cos(y)*1/f*sin(f*x+g*y)
> > > > > > > int(- cos((f-1)*x + g*y)  ./  4 ) = -1/(f-1)*sin((f-1)*x + g*y)/4
> > > > > > > int(- 2*cos(y).*(cos(y))=-2*x*cos^2(y)
> > > > > > > int(- 2*cos(y).*cos(x)) = -2*cos(y)*sin(x)
> > > > > > > Thus
> > > > > > > int(3 - 2*cos(y).*cos(f*x+g*y) - cos((f-1)*x + g*y)  ./  4 -
> > > > > > > 2*cos(y).*(cos(y) + cos(x)) dx) =
> > > > > > > 3*x-2*cos(y)*1/f*sin(f*x+g*y)-1/(f-1)*sin((f-1)*x + g*y)/
> > > > > > > 4-2*x*cos^2(y)-2*cos(y)*sin(x)
> > > > > > > Evaluate in the limits between -pi and pi:
> > > > > > > 3*pi-2*cos(y)*1/f*sin(f*pi+g*y)-1/(f-1)*sin((f-1)*pi + g*y)/
> > > > > > > 4-2*pi*cos^2(y) -
> > > > > > > ( 3*(-pi)-2*cos(y)*1/f*sin(-f*pi+g*y)-1/(f-1)*sin(-(f-1)*pi + g*y)/
> > > > > > > 4-2*(-pi)*cos^2(y)) =
> > > > > > > 6*pi-2*cos(y)*1/f*(sin(f*pi+g*y)-sin(-f*pi+g*y))-0.25/
> > > > > > > (f-1)*(sin((f-1)*pi + g*y)-sin(-(f-1)*pi + g*y))-4*pi*cos^2(y)
> > > > > > > Now use the
> > > > > > >http://integrals.wolfram.com
> > > > > > > to integrate this expression term by term.
> > > > > > > (Of course you have to substitute the y in the expressions by an x
> > > > > > > because the wolfram-integrator assumes
> > > > > > > the functions to integrate to depend on the variable x).
> > > > > > > Subsequently evaluate in the limits between -pi and pi.
> >
> > > > > > > The cases in which f=0 and |g|=1 may be treated seperately.
> >
> > > > > > > Best wishes
> > > > > > > Torsten.
> >
> > > > > > Im really sorry. I see where some of my misgivings have coome from.
> >
> > > > > > In my initial post I forgot to include the brackets around the denominator. The function is:
> > > > > > (3 - 2*cos(y).*cos(f*x+g*y) - cos((f-1)*x + g*y))  ./  (4 - 2*cos(y).*(cos(y) + cos(x))) dxdy
> >
> > > > > > Really sorry about that.- Zitierten Text ausblenden -
> >
> > > > > > - Zitierten Text anzeigen -
> >
> > > > > Then - as a further difficulty - you will have problems with the
> > > > > points where your
> > > > > denominator gets zero, e.g. (0,0) or (pi,pi) (where (pi,pi) is more
> > > > > critical than
> > > > > (0,0) because the numerator usually is different from 0 there).
> > > > > This might explain the error message of quad2d (if you entered the
> > > > > function
> > > > > correctly in the MATLAB-file).
> >
> > > > > Best wishes
> > > > > Torsten.
> >
> > > > Hi thanks for reply.
> >
> > > > Yes, this is the problem. Any ideas how I might somehow be able to manoeuvre this???
> >
> > > > cheers- Zitierten Text ausblenden -
> >
> > > > - Zitierten Text anzeigen -
> >
> > > Thinking about your function, I'm pretty sure that it is not
> > > integrable over
> > > [-pi;pi]x[-pi;pi] (at least for all values of f and g such that the
> > > numerator -
> > > evaluated at (+-pi/+-pi) - is different from 0).
> > > The reason is that in the neighbourhood of (+-pi/+-pi), the
> > > denominator
> > > behaves like (3*x^2+y^2) near (0,0), and the function 1/(3*x^2+y^2) is
> > > not
> > > integrable over a domain including (0,0).
> >
> > > Best wishes
> > > Torsten.
> >
> > Hi,
> >
> > Thanks for your reply and continual help - its much appreciated!
> >
> > -How can I get round this then?
> > -Could I integrate from -pi to (0-delta) and add the result to the integration from (0+delta) to pi??? Would that work?
> 
> I did not yet analyze if (0,0) is a problem - the problem I had in
> mind are the corner points of your domain
> (i.e. (pi/pi),(-pi/pi),(-pi,-pi),(pi,-pi)).
> It won't help to integrate up to a certain delta away from the
> singularity because the result will
> depend on that delta - simply because the integral of your function
> over [-pi;pi]x[-pi;pi] is infinity.
> 
> > -Additionally, I dont quite understand why matlab is happy to integrate the function for certain values of f and g and not others?
> >
> 
> I suspect these are numerical artefacts. Could you specify values of f
> and g for which the integration
> is successful and report the value MATLAB suggests for the integral ?
> 
> > merci- Zitierten Text ausblenden -
> >
> > - Zitierten Text anzeigen -
> 
> Best wishes
> Torsten.

Hi Thanks Torsten.

Ideally-------------
Ideally, I would like to evaluate the integral for all points where (f+g) is 

odd in the interval (0,0) to (2000,2000).

Test-------------
However, what I did:
Integrating from (0,0) to (0,1000) and (1,0) to (1,1000) and (2,0) to (2,000). 

Essentially taking the first 3 rows.
(quad2d(@(x,y)my_func(x,y,f,g),-pi,pi,-pi,pi,'MaxFunEvals',1e5));
And storing results in a 3x1001 array.
[*Note: ideally I want to run to 2000,2000 but will take along time. Also, I have included a check for (f+g) = even, and not computed that result.]
temp = (f+g);
if( floor((temp/2)) == (temp/2) )
flag_even=1;

Errors:-------------
Matlab outputs multiple errors of 2 kinds:
1. Warning: Reached the maximum number of function evaluations (100000). The result fails the global error test. 
> In quad2d at 248
  In Grid_Mesh_Hex at 94
2. Warning: Non-finite result. The integration was unsuccessful. Singularity likely. 
> In quad2d at 242
  In Grid_Mesh_Hex at 94

Result:-------------
For (f+g) is odd, all 3 rows give the correct answers up to and including column 133 for row 0. After that I get Inf up until 177 on row 0. Then more Inf until 241. Then a result. Then more Inf until 249 then a result. More Inf until 263.From 263 to 295 inclusive I get results again (results in the region of ~149).
This sort of sporadic output of Inf and numerical results continues. Similar sporadic patterns of results and Inf are evident for the other 2 rows. 

Example Results:
e.g) (1,0) = 26.31; (3,0) = 48.58; 
(0,1) = 26.31; (2,1) = 43.53
(1,2) = 48.58

Taking the curve fit of row 0 for the first 129 points (every other point where f+g is odd). The results follow the expected pattern of a*ln(x) + b.
For the higher values of f and/or g the change in the result should become less and less. These 'Inf' results are incorrect.
Hope this is useful.

Please help.

Thanks.