Double Integral Singularities quad2d

[S]

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!

[Torsten]

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.

[S]

Torsten <Torsten...@umsicht.fraunhofer.de> wrote in message <41a15ffe-e999-41b9...@l12g2000vby.googlegroups.com>...

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?

[Torsten]

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>...

> 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.

[S]

Torsten <Torsten...@umsicht.fraunhofer.de> wrote in message <db464627-33b0-456e...@f36g2000vbm.googlegroups.com>...

Thanks for your reply and so fast.

However, I do not have the symbolic toolbox installed :-(

Surely there must be another way?

[Torsten]

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>...

> 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.

[S]

Torsten <Torsten...@umsicht.fraunhofer.de> wrote in message <36eea4f4-4887-4557...@q13g2000vbd.googlegroups.com>...

Thanks for reply.

@ Torsten. Very good website btw. However, the website is unable to integrate the function.

[Roger Stafford]

"S" wrote in message <j838vr$9nu$1...@newscl01ah.mathworks.com>...

- - - - - - - - - - - -
Torsten is right. There is no point in subjecting matlab to a torturous numerical integration procedure for high values of f and g when a simple formula exists which can be obtained with pen and paper using elementary calculus. For example, the product 2*cos(y)*cos(f*x+g*y) can be manipulated this way:

2*cos(y)*cos(f*x+g*y) = 2*cos(y)*(cos(f*x)*cos(g*y)-sin(f*x)*sin(g*y)) =
cos(f*x)*2*cos(y)*cos(g*y) - sin(f*x)*2*cos(y)*sin(g*y) =
cos(f*x)*(cos((g+1)*y)+cos((g-1)*y)) ...
- sin(f*x)*(sin((g+1)*y)+sin((g-1)*y)) =
cos(f*x)*cos((g+1)*y) + cos(f*x)*cos((g-1)*y) ...
- sin(f*x)*sin((g+1)*y) - sin(f*x)*sin((g-1)*y)

In each of these last four terms the double integral can be obtained as the separate product of single integrals with respect to x and y respectively, and can therefore be easily evaluated.

Roger Stafford

[S]

"Roger Stafford" wrote in message <j840tv$s2i$1...@newscl01ah.mathworks.com>...

Thanks for your reply Roger.

Sorry, I dont follow you. Please forgive me for being a abit slow but I'm not sure what you are getting at.

[Roger Stafford]

"S" wrote in message <j842lp$4qi$1...@newscl01ah.mathworks.com>...

> Sorry, I dont follow you. Please forgive me for being a abit slow but I'm not sure what you are getting at.

- - - - - - - -

I'm not sure what puzzles you when you say "what you are getting at".

If it is the double integration analysis, consider for example the single term cos(f*x)*cos((g+1)*y) which I obtained. This expression has "separable" factors and its double integral from -pi to +pi for both x and y is therefore simply the product of the single integrals of cos(f*x) with x ranging from -pi to +pi and of cos((g+1)*y) with y ranging from -pi to +pi, which would be easily evaluated as:

(sin(pi*f)-sin(-pi*f))/f * (sin(pi*(g+1)-sin(-pi*(g+1))/(g+1) =
2*sin(pi*f)/f * 2*sin(pi*(g+1))/(g+1)

Thus its double integral has an analytic expression with no numerical integration required. The same reasoning applies for all the other terms that can be derived in your integrand.

If your uncertainty is due to not seeing why the numerical double integration would be "torturous", Torsten's suggestion that you should try plotting your function for high values of f and g should show you very dramatically why that would be a most difficult procedure. On an ordinary 'plot3' plot it would appear as a complete blur. It's no wonder that 'quad2d' refused to cooperate.

Roger Stafford

[Torsten]

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>...

> @ 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.

[S]

Torsten <Torsten...@umsicht.fraunhofer.de> wrote in message <272cde32-3068-4756...@c1g2000vbw.googlegroups.com>...

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.

[Torsten]

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>...

> 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.

[S]

Torsten <Torsten...@umsicht.fraunhofer.de> wrote in message <d15f286c-7e52-4d35...@p16g2000yqj.googlegroups.com>...

Hi thanks for reply.

Yes, this is the problem. Any ideas how I might somehow be able to manoeuvre this???

cheers

[Torsten]

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>...

> 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.

[S]

Torsten <Torsten...@umsicht.fraunhofer.de> wrote in message <05f56d39-0191-45b8...@r28g2000yqj.googlegroups.com>...

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?
-Additionally, I dont quite understand why matlab is happy to integrate the function for certain values of f and g and not others?

merci

[Torsten]

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>...

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.

[S]

Torsten <Torsten...@umsicht.fraunhofer.de> wrote in message <77e489a5-70a9-42a2...@eh5g2000vbb.googlegroups.com>...

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.

[Torsten]

On 26 Okt., 17:01, "S " <simaher2...@yahoo.co.uk> wrote:
> Torsten <Torsten.Hen...@umsicht.fraunhofer.de> wrote in message <77e489a5-70a9-42a2-a577-3db527134...@eh5g2000vbb.googlegroups.com>...

Why did you hide until now that you take f and g to be _integers_
which sum up to an _odd_ number ?
In case you would have provided this information straight away, it
would have been obvious that your
integral exists, and much of the discussion so far would have been
superfluous.
Anyway:
Because quad2d claims that there can be singularities on the boundary
of the domain of integration,
I suggest that you integrate your function over 4 domains seperately
and add the results, namely
[0,pi]x[0,pi], [-pi,0]x[0,pi], [0,pi]x[-pi,0] and [-pi,0]x[-pi,0].
This should at least prevent the second warning of a non-finite
result.
Choosing greater values for AbsTol and RelTol (especially for large
values of f and g)
should result in shorter computation times.

Best wishes
Torsten.

[S]

Torsten <Torsten...@umsicht.fraunhofer.de> wrote in message <a9b5312a-1c22-438c...@v8g2000vbe.googlegroups.com>...

My apologies Torsten. Initially, I thought I was just having issues with correctly using the matlab function and so didn't think that info was necessary :-/

Many thanks for your help - it is truly appreciated.

Could I just ask how you came to your conclusion. Just so that I can sort such problems out myself in the future.

And also because I have another function I want to integrate for (f+g) is even... f = (1-cos(f.*x+g.*y))./(4-2*cos(y).*(cos(y)+cos(x)));

Thanks again

[Torsten]

On 28 Okt., 10:03, "S " <simaher2...@yahoo.co.uk> wrote:
> Torsten <Torsten.Hen...@umsicht.fraunhofer.de> wrote in message <a9b5312a-1c22-438c-af29-b81f2238e...@v8g2000vbe.googlegroups.com>...

Which conclusion do you mean ?

Best wishes
Torsten.

[S]

Torsten <Torsten...@umsicht.fraunhofer.de> wrote in message <0dec8726-c843-401c...@j20g2000vby.googlegroups.com>...

I meant how you came to this:

[0,pi]x[0,pi], [-pi,0]x[0,pi], [0,pi]x[-pi,0] and [-pi,0]x[-pi,0].

I implemented it as follows:
** For f+g = Odd **
temp = (quad2d(@(x,y)my_func(x,y,f,g),0,pi,0,pi,'MaxFunEvals', 1e5));
temp = temp + (quad2d(@(x,y)my_func(x,y,f,g),-pi,0,0,pi,'MaxFunEvals', 1e5));
temp = temp + (quad2d(@(x,y)my_func(x,y,f,g),0,pi,-pi,0,'MaxFunEvals', 1e5));
temp = temp + (quad2d(@(x,y)my_func(x,y,f,g),-pi,0,-pi,0,'MaxFunEvals', 1e5));
result(f,g) = temp;

I ran for f=0 and f=1 and g up to 1000 and I still get errors of the sort:

Warning: Non-finite result. The integration was unsuccessful. Singularity likely.
> In quad2d at 242

In Grid_Mesh_Hex at 102

Warning: Non-finite result. The integration was unsuccessful. Singularity likely.
> In quad2d at 242

In Grid_Mesh_Hex at 103

Warning: Non-finite result. The integration was unsuccessful. Singularity likely.
> In quad2d at 242

In Grid_Mesh_Hex at 104

Warning: Non-finite result. The integration was unsuccessful. Singularity likely.
> In quad2d at 242

In Grid_Mesh_Hex at 105

[Torsten]

On 28 Okt., 15:54, "S " <simaher2...@yahoo.co.uk> wrote:
> Torsten <Torsten.Hen...@umsicht.fraunhofer.de> wrote in message <0dec8726-c843-401c-86cb-3026655dd...@j20g2000vby.googlegroups.com>...

I thought it would be a good idea to place the singularity in (0,0)
(which is in the interior of the domain)
on the boundary because the description of quad2d said that it could
handle singularities
there.
If it nevertheless refuses to integrate your function, I can't help
you.

If I were you, I'd just try the brute-force method:
Choose a grid as -pi=x1<x2<...<xn=pi, -pi=y(1)<y(2)<...<y(m)=pi in the
coordinate directions and
approximate your integral by
sum_i sum_j f((x(i)+x(i+1))/2,(y(j)+y(j+1))/2)*(x(i+1)-x(i))*(y(j+1)-
y(j)).
Choose the grid finer and finer until this sum converges to a stable
value.

Best wishes
Torsten.