Prev: How do I freeze the axes?
Next: ARMA estimated value prediction
From: Roger Stafford on 20 Mar 2010 14:10
"Vivek Saxena" <maverick280857(a)yahoo.com> wrote in message <ho24mk$mb1$1(a)fred.mathworks.com>...
> Hi Roger, I think I see your point. But what if I use quad2d?

If you use quad2d, then you could use the original x,y coordinates directly without any transformation, though that is a bit awkward because in general you have to break the integration into two separate ranges which depend on how the three points of the triangle are oriented. The triangle's edges may not be aligned with either of the coordinate axes.

It might be more convenient to use an affine transformation that transforms your triangle into another triangle that is oriented along the coordinate directions so that separate integration ranges would not be required. For example, you could say

x = x1 + u*(x2-x1) + v*(x3-x1)
y = y1 + u*(y2-y1) + v*(y3-y1)

This transforms the triangular region u >= 0, v >= 0, u+v <= 1 into your given x,y triangle. Thus your quad2d in the u,v domain could have simple limits: u from 0 to 1 and v from 0 to 1-u, and that is within quad2d's capabilities.

Of course you must still use the correct Jacobian factor in the integral if you do this transformation, but it is now just a constant quantity, namely the ratio between the two triangles' respective areas.

> So, if I understand you correctly, the triangle is mapped to the square {(u,v): 0<= u <= 1, 0 <= v <= 1} in (u,v) space.

Yes, that is correct (for the earlier transformation using dblquad.)

Roger Stafford
From: Vivek Saxena on 20 Mar 2010 23:02
>
> Yes, that is correct (for the earlier transformation using dblquad.)
>
> Roger Stafford

Hi again Roger. I thought about this a bit, and I'm unable to convince myself that the nonlinear transformation you wrote really does map a triangle into a square. I guess you're right, but when I did the integral, the values all came out wrong (I have an analytical solution to the final problem with which I can compare). And since this is the only additional term in my calculation, I have no reason to suspect the rest of the code which has been working so far.

Anyhow, I ran the following code

clear, close all;
clc;
x1 = 0, x2 = 1, x3 = 1, y1 = 0, y2 = 0, y3 = 1;
c = 1;
for u=0:0.01:1
for v=0:0.01:1
x(c) = (1-u)*x1 + u*((1-v)*x2+v*x3);
y(c) = (1-u)*y1 + u*((1-v)*y2+v*y3);
c = c + 1;
end
end
scatter(x,y);

That is I begin with a triangle which is right angled, with its vertices in (x, y) space as (0,0), (0,1) and (1,1) for simplicity. Now I use this transformation to construct x and y points for u and v ranging from 0 to 1. This gives me the coordinates of all the points inside the triangle (as required). But here, I have assumed that the corresponding region in (u, v) space is a square (because I allow u and v to run freely). So the forward transformation is a map from a square to a triangle. I am convinced about that. But how do I prove that the inverse transformation maps the triangle to the square {(u,v) : 0<=u<=1, 0<=v<=1}?

Perhaps this is a trivial question, but I can't seem to get my head around it. I am using dblquad, and my functions look like this:

function f = abcintegrand(u,v,a,b,c,k_0,x1,x2,x3,y1,y2,y3)

J = sqrt(-1);
argexp = (1-u).*x1 + u.*((1-v).*x2 + v.*x3);
jacob = (-x1+(1-v).*x2+v.*x3).*(-u.*y2+u.*y3)-(-y1+(1-v).*y2+v.*y3).*(-u.*x2+u.*x3);
p = (a + b*((1-u).*x1 + u.*((1-v).*x2 + v.*x3)) + c*((1-u).*y1 + u.*((1-v).*y2 + v.*y3))).*exp(-J*k_0*argexp);
f = p.*jacob;

FUNCTION CALL:

fun = @(u,v) abcintegrand(u,v,a(i),b(i),c(i),k_0,x1,x2,x3,y1,y2,y3);
b_integral(e,i) = (E_0*k_0*k_0*(eps_r(e)-1/mu_r(e))/(2*D(e)))*dblquad(fun,0,1,0,1);
From: Roger Stafford on 21 Mar 2010 20:39
"Vivek Saxena" <maverick280857(a)yahoo.com> wrote in message <ho427a$7d4$1(a)fred.mathworks.com>...

> ..... when I did the integral, the values all came out wrong .....

You will find that making a transformation of variables in evaluating multiple integrals using a Jacobian determinant is one of the fundamental methods in advanced calculus. It has been around a long time and you can surely count on it being valid. Note however that you must use the absolute value of the Jacobian determinant! I notice that you did not use the absolute value in your example. I should have made that clear earlier in this thread. See these three links:

http://en.wikipedia.org/wiki/Integration_by_substitution
http://wapedia.mobi/en/Integration_by_substitution
http://planetmath.org/encyclopedia/ChangeOfVariablesInIntegralOnMathbbRn.html

In your particular case of integration over a triangularly-shaped area, the transformation I gave you has a Jacobian determinant of:

u*det([x1,x2,x3;y1,y2,y3;1,1,1])

(Try simplifying your expression for 'jacob' and you will see.) However, as I said, what you need is the absolute value of the Jacobian as a factor, so if your three triangle vertices are oriented in clockwise order around the triangle, you would need to change the sign of your 'jacob' variable. You need its absolute value.

It is obvious that if your original integrand is unity, then its double integral over the triangle will of course be its area. The above Jacobian's absolute value is u multiplied by twice the triangle's area, so the integral in the u,v space will be int, u=1 to u=1, 2*area*u du = area also. The two integrals would thus be equal.

The same should be true of any integrand you might choose, including the one you used. I would suggest you try experimenting with very simple integrands to begin with, such as f(x,y) = x, y, x*y, or x^2+y^2 - something like that - to build up confidence in the method. With a little effort you should be able to find analytic expressions for both forms of the double integral in these simple cases and be able to show them to be identical.

If the problem is not due to a negative Jacobian factor, you should realize that dblquad works with a default error tolerance of only six decimal places unless you specify a smaller tolerance. How far off were your results?

> I'm unable to convince myself that the nonlinear transformation you wrote really does map a triangle into a square.

The inverse of the transformation

x = (1-u)*x1 + u*((1-v)*x2 + v*x3)
y = (1-u)*y1 + u*((1-v)*y2 + v*y3)

is this (courtesy of matlab symbolic toolbox)

t = ((x-x1)*(y3-y2)-(y-y1)*(x3-x2)); % Temporary variable
u = t/det([x1,x2,x3;y1,y2,y3;1,1,1]);
v = (y*(x2-x1)-x*(y2-y1)+x1*y2-x2*y1)/t;

Therefore the following code should demonstrate that this inverse does indeed transform all of the interior of a triangle in x,y into the interior of the unit cube in u,v. I have taken the liberty of using my FEX program, randfixedsum, to obtain barycentric coordinates for creating points within a random triangle. It is available at

http://www.mathworks.com/matlabcentral/fileexchange/9700

Here is the code:

x1 = rand-1.1; y1 = rand; % Choose a random triangle off to the left
x2 = rand-1.1; y2 = rand;
x3 = rand-1.1; y3 = rand;
r = randfixedsum(3,2000,1,0,1); % All 2000 columns sum to 1
a = r(1,:); b = r(2,:); c = r(3,:);
x = a*x1+b*x2+c*x3; % a,b,c are barycentric coord.
y = a*y1+b*y2+c*y3;
t = ((x-x1)*(y3-y2)-(y-y1)*(x3-x2)); % Here is the inverse transform
u = t/det([x1,x2,x3;y1,y2,y3;1,1,1]);
v = (y*(x2-x1)-x*(y2-y1)+x1*y2-x2*y1)./t;
plot(x,y,'y.',[x1,x2,x3,x1],[y1,y2,y3,y1],'b-',...
u,v,'r.',[0,1,1,0,0],[0,0,1,1,0],'b-')
axis equal

As you can see, it transforms uniformly distributed points within the triangle over to points within the unit square in u,v, though no longer uniformly spaced. This non-uniformity reflects the variation of the Jacobian factor - the density is smaller as u becomes smaller on the left side.

Roger Stafford
From: Vivek Saxena on 21 Mar 2010 21:14
Hi Roger, thanks for the detailed reply. I am familiar with the theory of the method (being an engineering student myself :-) ). That is not the problem. I was wondering however, that given the forward transformation (x, y) --> (u, v) which clearly maps the triangle into the square (simply because u and v are now independent), can we show that (u,v) --> (x,y) is unique. Anyway since my last post, I managed to convince myself about this. That was just my way of trying to motivate the choice of the transformation.

Also, when I said that the values of the integral being computed are turning out to be wrong, I was referring not to the transformation but rather to my usage of dblquad. I'm inclined to believe that something is wrong with my usage of that function. I should've been more explicit. (Thanks for your help btw.)
From: Bruno Luong on 22 Mar 2010 03:20
"Vivek Saxena" <maverick280857(a)yahoo.com> wrote in message <ho6g8s$os1$1(a)fred.mathworks.com>...
> Hi Roger, thanks for the detailed reply. I am familiar with the theory of the method (being an engineering student myself :-) ). That is not the problem. I was wondering however, that given the forward transformation (x, y) --> (u, v) which clearly maps the triangle into the square (simply because u and v are now independent), can we show that (u,v) --> (x,y) is unique. Anyway since my last post, I managed to convince myself about this. That was just my way of trying to motivate the choice of the transformation.
>
> Also, when I said that the values of the integral being computed are turning out to be wrong, I was referring not to the transformation but rather to my usage of dblquad. I'm inclined to believe that something is wrong with my usage of that function. I should've been more explicit. (Thanks for your help btw.)

One way to check the correctness of the parametrization is to integrate a polynomial P(x) on the triangle. If P(x) the polynomial of
- order 0: the integral result is P(x1)*|T|
- order 1; the integral result is (P(x1)+P(x2)+P(x3)) * |T|/3
etc...

x1, x2, x3 are three corners, and |T| is the area of triangle

Once the implementation gives correct result, you can then replace P(x) by your function.

Bruno
First  |  Prev  |  Next  |  Last
Pages: 1 2 3 4
Prev: How do I freeze the axes?
Next: ARMA estimated value prediction