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