COMPGROUPS.NET | Search | Post Question | Groups | Stream | About | Register

Help with discrete double integral

• Follow

```Hi everyone,

My problem is : I want to integrate f(x,y)*q(x) over x&#65374;(a,b),y~(c,d).
f(x,y) and q(x)  are discrete numbers not functions.
The f(x,y) is a two-dimensional function which is dispalyed in discrete numbers. I couldn't figure out how to do this integration. Any suggestion is appreciated.

Chen
```
 0

```I'm not sure I understand this.  How are these numbers generated?  You need to have some knowledge of the problem before applying numerical methods to it.

"Chen Chen" <chenchen_ee@yahoo.com> wrote in message <i3p629\$8oh\$1@fred.mathworks.com>...
> Hi everyone,
>
>    My problem is : I want to integrate f(x,y)*q(x) over x&#65374;(a,b),y~(c,d).
>    f(x,y) and q(x)  are discrete numbers not functions.
>    The f(x,y) is a two-dimensional function which is dispalyed in discrete numbers. I couldn't figure out how to do this integration. Any suggestion is appreciated.
>
> Chen
```
 0

```"Joseph " <don'twannapostit@nopers.com> wrote in message <i3p9ap\$9ih\$1@fred.mathworks.com>...
> I'm not sure I understand this.  How are these numbers generated?  You need to have some knowledge of the problem before applying numerical methods to it.

Hi Joseph,

I just simplify the problem i want to solve to this imple double mathematical problem. I thought it would be confusing if I explained the whole thing.
Sorry I didn't state the problem clearly. Here is a specific version.
These numbers are given sets of data. That is: x and y are two sets of data which are given at the first place. Let's say both the x and y are 1 X N vectors. f(x,y) is given too and it's an N X N matrix. q(x) is a function depending on x only and I know the function already. So q(x) is a 1 X N vectors with known numbers.
My problem is that I don't know how to integrate f(x,y)*f(x) over x and y numerically. I hope I expained it clearly this time.

Chen
```
 0

```Chen Chen wrote:
> "Joseph " <don'twannapostit@nopers.com> wrote in message
> <i3p9ap\$9ih\$1@fred.mathworks.com>...
>> I'm not sure I understand this.  How are these numbers generated?  You
>> need to have some knowledge of the problem before applying numerical
>> methods to it.
>
>
> Hi Joseph,
>
>    I just simplify the problem i want to solve to this imple double
> mathematical problem. I thought it would be confusing if I explained the
> whole thing.    Sorry I didn't state the problem clearly. Here is a
> specific version.
>    These numbers are given sets of data. That is: x and y are two sets
> of data which are given at the first place. Let's say both the x and y
> are 1 X N vectors. f(x,y) is given too and it's an N X N matrix. q(x) is
> a function depending on x only and I know the function already. So q(x)
> is a 1 X N vectors with known numbers.    My problem is that I don't
> know how to integrate f(x,y)*f(x) over x and y numerically. I hope I
> expained it clearly this time.

Typical numeric integrals rely upon implied continuity, that f(x) is some
value then there is some (unknown but present) path from f(x-1) to f(x) to
f(x+1). When, though, you say that the data is discrete, then you disclaim any
implied continuity, leaving us with no knowledge of what you mean by
"integrate" in this situation. Do you want each value f(x,y) to apply over the
semi-closed intervals [x, x+1) and [y, y+1), or do you want f(x,y) to apply
over [x-1/2, x+1/2) and (y-1/2, y+1/2) or something else? And in the integral,
should the neighbouring values be ignored when considering f(x,y)*q(x), or
should the neighbours be considered?

With the information you have given, it sounds to me as if you might simply want

sum( reshape( bsxfun('times', q(:), f), [], 1) )
```
 0

```Hi walter,

The x and y are given domain. Let's say x~(-10,10), y ~(10,100).

Chen
```
 0

```Chen wrote:
>      The x and y are given domain. Let's say x~(-10,10), y ~(10,100).

That's nice, but it doesn't answer even one of the questions I asked.

Does the N x N matrix of values for f(x,y) *define* f(x,y), or does the matrix
represent a *sampling* of a function at discrete grid points? If it represents
a *sampling* of a function, then is the function known to have continuous
first derivatives? Is it known to have continuous second derivatives as well?

Your matrix for f(x,y) is N x N, but the x and y ranges you give are not the
same: one has a range of 20 and the other has a range of 90. That implies
finding the area of rectangles rather than squares, a fact which would affect
numerical integration formula. You are using the ~ operator, which in
statistics means "is distributed according to", but you have a bunch of
(x,y,f) triples as implied by "is distributed according to" then you would not
be able to put them into a regular NxN grid (well, not meaningfully). But
still you used the ~ statistical operator: is the implication that you _do_
have a rectangular grid but that the grid spacing is irregular?
```
 0

```Hi walter,

Your questions are really good. Those conditions should be taken care of when you are doing puring mathematical calculation (1st and 2nd derivative). But I am only doing this for engineering use. The functions are basicly well behaved .
Sorry about the wrong use of ~. when I said x~(-10,10), I meant that x is from -10 to 10.
Let's put it this way. I want to integrate a function f(x,y)*p(x) over x and y.
(1)
If f(x,y) = cos(x)*sin(y) and p(x) = cos(2*x) (!! they are functions here ); .    Theoretically, I can integrate this with no trouble.
(2)
But now the thing is that I need to integrate discrete data set of f(x,y) and p(x) (!! They are discrete numbers here). and both x and y are 1 X N vectors. Then the problem becomes tricky to me. Now the f(x,y) is a N X N matrix and p(x) is a 1 X N vector.

I hope you understand what i am saying..
```
 0

```"Chen Chen" <chenchen_ee@yahoo.com> wrote in message <i3pfrd\$aca\$1@fred.mathworks.com>...
>     I just simplify the problem i want to solve to this imple double mathematical problem. I thought it would be confusing if I explained the whole thing.
>     Sorry I didn't state the problem clearly. Here is a specific version.
>     These numbers are given sets of data. That is: x and y are two sets of data which are given at the first place. Let's say both the x and y are 1 X N vectors. f(x,y) is given too and it's an N X N matrix. q(x) is a function depending on x only and I know the function already. So q(x) is a 1 X N vectors with known numbers.
>     My problem is that I don't know how to integrate f(x,y)*f(x) over x and y numerically. I hope I expained it clearly this time.
>
> Chen
- - - - - - - - - -
With discrete-valued variables you cannot use any of matlab's "quad" type integration routines which require functions.  However you can use an iterated form of trapz which does trapezoidal integration.  Let X be a column vector of monotone x values in which X(1) = a and X(N) = b.  Let Y be a column vector of monotone y values for which Y(1) = c and Y(M)= d.  Let F be an N x M array in which

F(i,j) = f(X(i),Y(j)).  (Your f(x,y) )

Let Q be a column vector in which Q(i) = q(X(i)) (your q(x) ).  Then do this:

I = trapz(X,Q.*trapz(Y,F.').');

This will give the two-dimensional "trapezoidal" approximation to your desired double integral over the stated rectangle.

The approximation made is this.  In each rectangular cell with corners at

(X(i),Y(j), (X(i+1),Y(j), (X(i),Y(j+1), and (X(i+1),Y(j+1),

the integral is approximated by the average of the integrand q(x)*f(x,y) evaluated at these four corners multiplied by the (signed) area of the rectangle.  The final value I is the sum of these.

This integration can be considered as a first-order approximation to the exact double integral.  I think if you look hard enough there may be some higher order integration routines listed on the file exchange which could possess higher accuracy, though if so, I don't remember where they are located.

Roger Stafford
```
 0

```Chen wrote:
> Hi walter,
>
>    Your questions are really good. Those conditions should be taken care
> of when you are doing puring mathematical calculation (1st and 2nd
> derivative). But I am only doing this for engineering use.

Which numeric integration routines you can use is affected by how many layers
of continuity you assume.
```
 0

```Hi Roger,

I tried this:
f(x,y) = cos(x).*sin(y) and q(x) = cos(x).x is from 0 to pi/2,y is from pi/2 to 2*pi. The exact integral value of f(x,y)*q(x) over x and y is - pi/4.

The code is :
%%% Integration test

x = linspace(0,pi/2,100);
y = linspace(pi/2,2*pi,100);

func = @(x,y) cos(x).^2.*sin(y);

% method 2,using double trapz
[X Y] = meshgrid(x,y);
f = cos(X).*sin(Y);
q = cos(x);
I2 = trapz(x,q.*trapz(y,f.').');

%%%%% End of code %%%%%%%%%

The results are: I1 = -0.7854(which is - pi/4) and I2 = 5.2042e-017 (which is zero actually).

Did I understand your method wrong?
Thanks a lot

-Chen
```
 0

```"Chen " <neversaynever@never.org> wrote in message <i3qfon\$2h5\$1@fred.mathworks.com>...
> Hi Roger,
>
>    Thank you for your reply. I tried your method and it seems like not working right.
> I tried this:
> f(x,y) = cos(x).*sin(y) and q(x) = cos(x).x is from 0 to pi/2,y is from pi/2 to 2*pi. The exact integral value of f(x,y)*q(x) over x and y is - pi/4.
>
>   The code is :
> %%% Integration test
>
> x = linspace(0,pi/2,100);
> y = linspace(pi/2,2*pi,100);
>
> func = @(x,y) cos(x).^2.*sin(y);
>
> % method 2,using double trapz
> [X Y] = meshgrid(x,y);
> f = cos(X).*sin(Y);
> q = cos(x);
> I2 = trapz(x,q.*trapz(y,f.').');
>
> %%%%% End of code %%%%%%%%%
>
> The results are: I1 = -0.7854(which is - pi/4) and I2 = 5.2042e-017 (which is zero actually).
>
> Did I understand your method wrong?
> Thanks a lot
>
> -Chen
- - - - - - - - - - -
The problem is in your use of meshgrid.  If you recall, my definition of F was that  F(i,j) = f(X(i),Y(j)).  However when you say

x = linspace(0,pi/2);
y = linspace(pi/2,2*pi);
[X Y] = meshgrid(x,y);
f = cos(X).*sin(Y);

that makes f backwards so that f(i,j) = cos(x(j))*sin(y(i)).  In other words your f is the transpose of the one I defined.

If you write

x = linspace(0,pi/2).';
y = linspace(pi/2,2*pi).';
[Y,X] = meshgrid(y,x);
f = cos(X).*sin(Y);
q = cos(x);
I = trapz(x,q.*trapz(y,f.').');

it should work properly.

Or if you wish you can write

x = linspace(0,pi/2);
y = linspace(pi/2,2*pi);
[X,Y] = meshgrid(x,y);
f = cos(X).*sin(Y);
q = cos(x);
I = trapz(x,q.*trapz(y,f));

and that should also work.

Of course you must realize that this double trapz scheme cannot achieve the same accuracy as quad2d for explicitly defined functions since the latter has unlimited access to such functions, unless you use a sufficiently large number of points in your trapz mesh arrays.

Roger
```
 0
Reply ellieandrogerxyzzy (4732) 8/10/2010 5:46:04 AM

```Hi Roger,

I see where the problem is. Thanks a lot !

-Chen
```
 0

```Hi Roger,

I still have one problem to bother you.
Based on the thing you explained to me yesterday, what if i want to integrat
f(x,y).*p(x).*q(y) over x and y? where f(x,y) is still a N X N matrix, p(x) and q(y) are 1 X N vectors.

-Chen
```
 0

```"Chen " <neversaynever@never.org> wrote in message <i3s7a9\$bdo\$1@fred.mathworks.com>...
> Hi Roger,
>
>    I still have one problem to bother you.
>    Based on the thing you explained to me yesterday, what if i want to integrat
>    f(x,y).*p(x).*q(y) over x and y? where f(x,y) is still a N X N matrix, p(x) and q(y) are 1 X N vectors.
>
> -Chen
- - - - - - - - - - -
I'll put it in words.  Suppose there are n elements in x and m elements in y.  If you want to integrate with respect to y in the inner iterated integral, you need an m by n array in which the element in the k-th row amd h-th column is equal to q(y(k))*f(x(h),y(k)).  (You can get that using bsxfun or repmat.)  That would be the second argument in the inner trapz.  The first argument should be a vector whose k-th element is y(k).  The result from trapz would be an n-element row vector in which the h-th element would be the (approximate) integral of q(y)*f(x,y) with respect to y taken along the h-th column at the fixed value x = x(h).

In the outer integral the second argument should be an n-element vector consisting of each of these n result elements multiplied by the corresponding p(x(h)).  The first argument should be a vector in which the h-th element is x(h).  The result of this would be the (approximate) double integral over the corresponding rectangle in x and y.

I will leave it to you to write the necessary matlab code to accomplish all this.  My advice is to first try it on an f, p, and q for which you know the exact double integral, as you did before, and compare it to the trapz result.  You should try it with m and n set at different values from one another.  Depending on how large m and n are, you should get a reasonably accurate matchup.  If not, you should check for errors in the coding.

Also before you commence coding I would advise a careful reading of the help file for trapz, especially for the case where the second argument is a matrix.  Try it out for a few very simple cases with very few points to get a feeling for how it works.  You should also try verifying with for-loops the claim I made earlier that the result of the double trapz would be the sum over all the individual rectangle cells of the area of the cell times the average of the integrand at its four corners.  The only differences here should be those due to round-off errors.

Note: You can also do this double integration by incorporating the entire integrand in the m by n matrix for the second argument of the inner trapz.  Its k-th row amd h-th column should be equal to p(x(h))*f(x(h),y(k))*q(y(k)), and again bsxfun or repmat would be handy in creating it.  The second argument of the outer trapz would then be just the result from the inner trapz.

Roger Stafford
```
 0

```Hi Roger,

You are such a great teacher. I tried your method and it worked!

x = linspace(0,pi/2);
y = linspace(pi/2,2*pi);

func = @(x,y) sin(x).*cos(y).*cos(x).*sin(y);

% method 2,using double trapz
[X Y] = meshgrid(x,y);
f = sin(X).*cos(Y);
f = f';
p = cos(x);
q = sin(y);
for i = 1: length(x)
for j = 1:length(y)
ff(i,j) = p(j)*f(j,i)*q(i);
end
end;
I2 = trapz(y,trapz(x,ff))
%%%%%%%%%%%%%%%%%%% end of test code %%%%%%%%%%5

if 100 points are tried, the results are I1 = - 0.2500; I2 = -0.2498
if 1000 points are tried,                    I1 = - 0.2500; I2 = -0.2500

Thank you very much and I appreciate it.

-Chen
```
 0

```i don't know how to use double integration for my formula  which is

A(m,n)= double integration of  f(x,y)*cos(m*theta)*cos(2*n*pi*x/L) dy dx in the range of ymin=0 ymax=2*pi*r xmax=L xmin=0

where f(x,y) is a discrete datas in the range x=16 and y=73
theta=0:5:2pi and n=1 to 17 and m=1 to 73 , L=150mm r-47.5mm

can anybody help me how to solve this problem..

0
```
 0

```"prabhakaran m" <prabha.gahon@gmail.com> wrote in message <i726mh\$k3h\$1@fred.mathworks.com>...
> i don't know how to use double integration for my formula  which is
>
> A(m,n)= double integration of  f(x,y)*cos(m*theta)*cos(2*n*pi*x/L) dy dx in the range of ymin=0 ymax=2*pi*r xmax=L xmin=0
>
> where f(x,y) is a discrete datas in the range x=16 and y=73
> theta=0:5:2pi and n=1 to 17 and m=1 to 73 , L=150mm r-47.5mm
>
> can anybody help me how to solve this problem..
- - - - - - - -
Prabhakaran, I refer you to the eighth (and fourteenth) articles in this thread in which I explain how to doubly integrate this type of discrete valued integrand using trapz.  Just read those articles and then read the documentation for trapz very carefully.  There is no point in me (or anyone else) explaining the same thing all over again.

Roger Stafford
```
 0

16 Replies
548 Views

Similiar Articles:

7/23/2012 1:21:13 AM