draw an Ellipsoid given its quaderatic equation in Matlab

Question

How can I draw an ellipsoid in Matlab with following equation:

$aX^2 + bY^2 + cZ^2 + d XY+e XZ +f YZ +gX +hY+iZ=R^2$

with (for example) $a=1; b=2; c=3; d=4; e=5; f=6; g=7; h=8; i=9; R^2=10.$

Answer 1

The equation you give is of the form $\alpha Z^2 + \beta Z + \gamma = 0$, so given the values of all the coefficients along with grid of $(X,Y)$ pairs, we can use the quadratic formula to solve for the two values of $Z$ that satisfy the equation:

[x,y] = meshgrid( (-50:0.1:50), (-50:0.1:50)' );

quadsoln1 = @(a,b,c)( ( -b + sqrt(b.^2 - 4*a.*c) ) ./ (2*a) );
quadsoln2 = @(a,b,c)( ( -b - sqrt(b.^2 - 4*a.*c) ) ./ (2*a) );

a = 1;
b = 2;
c = 3;
d = 4;
e = 5;
f = 6;
g = 7;
h = 8;
i = 9;
R = sqrt(10);

z1 = quadsoln1( c, e*x + f*y + i, a*x.^2 + b*y.^2 + d*x.*y + g*x + h*y - R );
z2 = quadsoln2( c, e*x + f*y + i, a*x.^2 + b*y.^2 + d*x.*y + g*x + h*y - R );

To plot the above, you can use plot3 to plot points or try something like mesh or surf. The latter two, however, are quite memory intensive and seem just unstable in general.

Answer 2

You can treat the ellipsoid as a transformed sphere, with points $v$ given by an equation of the form $v^t M v = 1$. Because the matrix $M$ will turn out to be symmetric with nonnegative eigenvalues, it has a square root, say $S$, so that $S^t S = M$. Then the equation becomes $v^t S^t S v = 1$, so $Sv$ is a vector $u$ satisfying $u^t u = 1$, i.e., a point of the unit sphere. Letting $v = S^{-1}u$, you're done. here's the code:

function ellipsoid(a, b, c, d, e, f, g, h, j, R)
M = (1/R^2) * [
    a d e g;
    0 b f h;
    0 0 c j;
    0 0 0 1];
M = (M + M')/2;
S = M^0.5

n = 30;
[th, ph] = meshgrid(linspace(0, 2*pi, n), linspace(0, pi, n));
spherepoints = zeros(4, n, n);
spherepoints(1, :,:) = cos(th) .* sin(ph);
spherepoints(2, :,:) = sin(th) .* sin(ph);
spherepoints(3, :,:) = cos(ph);
spherepoints(4, :, :) = 1;
sp = reshape(spherepoints, 4, []); 
epoints = S \ sp;
epoints = reshape(epoints, 4, n, n);

surf(squeeze(epoints(1, :,:)),...
    squeeze(epoints(2,:,:)),...
    squeeze(epoints(3,:,:)));
set(gca, 'DataAspectRatio', [1,1,1]);
figure(gcf);

Thanks to the anonymous reviewer who noted that I'd written epoints = S * sp instead of S \ sp.