Prev: Transparency in 2-D Area Plots
Next: Simulating an inverter
From: Roger Stafford on 10 Aug 2010 01:46
"Chen " <neversaynever(a)never.org> wrote in message <i3qfon$2h5$1(a)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);
>
> % method 1,using quad2d
> func = @(x,y) cos(x).^2.*sin(y);
> I1 = quad2d(func,0,pi/2,pi/2,2*pi);
>
> % 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
From: Chen on 10 Aug 2010 08:40
Hi Roger,

I see where the problem is. Thanks a lot !

-Chen
From: Chen on 10 Aug 2010 14:55
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
From: Roger Stafford on 10 Aug 2010 16:27
"Chen " <neversaynever(a)never.org> wrote in message <i3s7a9$bdo$1(a)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
From: Chen on 10 Aug 2010 23:14
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);

% method 1,using quad2d
func = @(x,y) sin(x).*cos(y).*cos(x).*sin(y);
I1 = quad2d(func,0,pi/2,pi/2,2*pi)

% 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
First  |  Prev  | 
Pages: 1 2 3
Prev: Transparency in 2-D Area Plots
Next: Simulating an inverter