Computing circumcenter of triangle in 2D with MATLAB I'm writing a finite volume program over a 2D triangular mesh, and at one point I need to calculate the circumcenters of the triangles. The equation given in class and that on Wikipedia give different results.
In class:
$$ \vec{r}_c = \vec{r}_i + \frac{( {\lVert \vec{r}_{ik} \rVert}^2 \vec{r}_{ij} - {\lVert \vec{r}_{ij} \rVert}^2 \vec{r}_{ik}) \times \vec{r}_{ik} \times \vec{r}_{ij}} { 2 {\lVert \vec{r}_{ik} \times \vec{r}_{ij} \rVert}^2} $$
Wikipedia:
$$ x_c = \frac{{\lVert A \rVert}^2(B_y-C_y)+{\lVert B \rVert}^2(C_y-A_y)+{\lVert C \rVert}^2(A_y-B_y)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}$$
$$y_c = \frac{{\lVert A \rVert}^2(C_x-B_x)+{\lVert B \rVert}^2(A_x-C_x)+{\lVert C \rVert}^2(B_x-A_x)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))} $$
Where $\vec{r}_{ab} = \vec{r}_b - \vec{r}_a$, and $(\vec{r}_i,\vec{r}_j,\vec{r}_k)$ and $(A,B,C)$ are the triangle vertices.
I implented both after running into problems with the final solution (cross section of fluid velocity in a strangely-shaped tube), and found they gave different solutions.
Plotting the triangles and calculated circumcenters gives the following.
Using in-class formula: 
Using Wikipedia formula: 
Circumcenters were calculated with the following lines of code.
In-class formula:
C = [I,0]+(cross(cross(((norm(K-I))^2*[J-I,0]-(norm(J-I))^2*[K-I,0]),[K-I,0]),[J-I,0]))/(2*(norm(cross([K-I,0],[J-I,0])))^2);
C = C(1,1:2);

Wikipedia:
D = 2*(I(1)*(J(2)-K(2))+J(1)*(K(2)-I(2))+K(1)*(I(2)-J(2)));
xc = ((I(2)^2+I(1)^2)*(J(2)-K(2))+(J(2)^2+J(1)^2)*(K(2)-I(2))+(K(2)^2+K(1)^2)*(I(2)-J(2)))/D;
yc = ((I(1)^2+I(2)^2)*(J(1)-K(1))+(J(1)^2+J(2)^2)*(K(1)-I(1))+(K(1)^2+K(2)^2)*(I(1)-J(1)))/D;
C = [xc -yc];

I figure that there must be something wrong with my implementation of the in-class formula, as many acute triangles don't have a calculated circumcenter inside of them, but I can't find the error. Can you?
PS: Yeah, I see the more elegant methods of implementing the Wikipedia formula. Meh.
 A: Substituting variables $(\vec{r}_i,\vec{r}_j,\vec{r}_k) = (A,B,C)$, with all z-coordinates equal to zero, into the 'in-class' formula gives the following:
$$\vec{r}_c = A+\frac{({\lVert C-A \rVert}^2(B-A)-{\lVert B-A \rVert}^2(C-A)) \times (C-A) \times (B-A)}{2{\lVert (C-A) \times (B-A) \rVert}^2}$$
$$\vec{r}_c = A+\frac{(({\lVert A \rVert}^2+{\lVert C \rVert}^2-2A \cdot C)(B-A)-({\lVert A \rVert}^2+{\lVert B \rVert}^2-2A \cdot B)(C-A)) \times (C-A) \times (B-A)}{2{\lVert (C-A) \times (B-A) \rVert}^2}  $$
$$\vec{r}_c = A+\frac{\left({\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A) -2\left[(A\cdot C)(B-A)-(A\cdot B)(C-A)\right]\right) \times (C-A) \times (B-A)}{2{\lVert (C-A) \times (B-A) \rVert}^2} $$
Taking a hint from the denominator, I assume the cross product follows the following order. (The terms to the right combine to form a unit vector.)
$$\vec{r}_c = A+\frac{{\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A) -2\left[(A\cdot C)(B-A)-(A\cdot B)(C-A)\right]}{2{\lVert (C-A)\times(B-A) \rVert}} \times \frac{(C-A)\times(B-A)}{\lVert (C-A)\times(B-A)\rVert}$$
$$\vec{r}_c = A+\left(\frac{{\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}-\frac{(A\cdot C)(B-A)-(A\cdot B)(C-A)}{\lVert (C-A)\times(B-A)\rVert}\right)\times \hat z$$
Let $Y=B-A$, $Z=C-A$, $A\cdot B=X\cdot Y$, and $A\cdot C=X\cdot Z$, giving the following:
$$\vec{r}_c = A+\left(\frac{{\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}-\frac{(X\cdot Z)Y-(X\cdot Y)Z}{\lVert (C-A)\times(B-A)\rVert}\right)\times \hat z$$
$$\vec{r}_c = A+\left(\frac{{\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}-\frac{X\times(Y\times Z)}{\lVert (C-A)\times(B-A)\rVert}\right)\times \hat z$$
$$\vec{r}_c = A+\left(\frac{{\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}+X\times\frac{(C-A)\times(B-A)}{\lVert (C-A)\times(B-A)\rVert}\right)\times \hat z$$
$$\vec{r}_c = A+(X\times\hat z)\times\hat z +\frac{{\lVert A \rVert}^2 (C-B) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}\times\hat z$$
Now, if I knew how to solve for $-X$, I'd show the following:
$$\vec{r}_c = \frac{{\lVert A \rVert}^2 (C-B) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}\times\hat z$$
In other words, the two formulas are off by a factor of -1, and it appears to be the Wikipedia one that maps improperly, as I had been 'fixing' it by setting yc=-yc -- the x-axis is symmetric, so I didn't notice that it was also flipped.
The 'in-class' formula needs to be implemented like this: (@Chris Taylor hit it on the head!)
C = [I,0]+(cross((norm(K-I))^2*[J-I,0]-(norm(J-I))^2*[K-I,0],cross([K-I,0],[J-I,0])))/(2*(norm(cross([K-I,0],[J-I,0])))^2);
C = C(1,1:2);

I'll have to try it out later, as it's time to head to class!
