1

Could you please help me with the following question: I want to solve a second order equation with two unknowns and use the results to plot an ellipse. Here is my function:

fun = @(x) [x(1) x(2)]*V*[x(1) x(2)]'-c

V is 2x2 symmetric matrix, c is a positive constant and there are two unknowns, x1 and x2. If I solve the equation using fsolve, I notice that the solution is very sensitive to the initial values

fsolve(fun, [1 1])

Is it possible to get the solution to this equation without providing an exact starting value, but rather a range? For example, I would like to see the possible combinations for x1, x2 \in (-4,4)

Using ezplot I obtain the desired graphical output, but not the solution of the equation. 

fh= @(x1,x2) [x1 x2]*V*[x1 x2]'-c;
ezplot(fh)
axis equal

Is there a way to have both? Thanks a lot!

Astrid
  • 43
  • 1
  • 2
  • 7

1 Answers1

3

you can take the XData and YData from ezplot:

c = rand;
V = rand(2);
V = V + V';
fh= @(x1,x2) [x1 x2]*V*[x1 x2]'-c;
h = ezplot(fh,[-4,4,-4,4]); % plot in range
axis equal
fun = @(x) [x(1) x(2)]*V*[x(1) x(2)]'-c;
X = fsolve(fun, [1 1]); % specific solution
hold on;
plot(x(1),x(2),'or');
% possible solutions in range
x1 = h.XData;
x2 = h.YData;

or you can use vector input to fsolve:

c = rand;
V = rand(2);
V = V + V';
x1 = linspace(-4,4,100)';
fun2 = @(x2) sum(([x1 x2]*V).*[x1 x2],2)-c;
x2 = fsolve(fun2, ones(size(x1))); 
% remove invalid values
tol = 1e-2;
x2(abs(fun2(x2)) > tol) = nan;

plot(x1,x2,'.b')

However, the easiest and most straight forward approach is to rearrange the ellipse matrix form in a quadratic equation form:

k = rand;
V = rand(2);
V = V + V';
a = V(1,1);
b = V(1,2);
c = V(2,2);
% rearange terms in the form of quadratic equation:
%     a*x1^2 + (2*b*x2)*x1 + (c*x2^2) = k;
%     a*x1^2 + (2*b*x2)*x1 + (c*x2^2 - k) = 0;
x2 = linspace(-4,4,1000);
A = a;
B = (2*b*x2);
C = (c*x2.^2 - k);
% solve regular quadratic equation
dicriminant = B.^2 - 4*A.*C;
x1_1 = (-B - sqrt(dicriminant))./(2*A);
x1_2 = (-B + sqrt(dicriminant))./(2*A);
x1_1(dicriminant < 0) = nan;
x1_2(dicriminant < 0) = nan;
% plot
plot(x1_1,x2,'.b')
hold on
plot(x1_2,x2,'.g')
hold off
user2999345
  • 4,195
  • 1
  • 13
  • 20
  • Thank you very much for your answer! In the first solution, I can retrieve the XData and YData, but how do I know for which pairs the solution is zero? In the second solution, the plot does not look like an ellipse, any idea how can I retireve an ellipse? – Astrid Apr 06 '17 at 14:02
  • 1
    `ezplot` plots **only** for `x`.`y` pairs which their solution is zero, so - **all of them**. because of the random matrices I used in this example you sometimes get an ellipse as a solution and sometimes asymptotic lines. use your own matrix and you should get an ellipse (in case it is the right solution). – user2999345 Apr 06 '17 at 14:11
  • The XData and YData both go from -4 to 4 in some step. From the elipse, it can be seen that e.g. -3.96 and -3.96 it's not a viable solution. So I would like to know which is the y corresponding to x = -3.96. – Astrid Apr 06 '17 at 14:16
  • The second method is very useful, but it is highly sensitive to the initial values, i.e. if instead of ones(size(x1)) I put zeros(size(x1)), the solution is very different. – Astrid Apr 06 '17 at 14:17
  • 1
    note the each `x1` value has two `x2` values which satisfy the equation (symmetric w.r.t some line), so you might just get one of the each time you run the code. I edited my second part of the answer to remove invalid values (which do not satisfy the equation). – user2999345 Apr 06 '17 at 14:24
  • Exactly, you are right, for each x1, I will have two possible x2. So there is no way I can retrieve the pairs from the ellipse solution? – Astrid Apr 06 '17 at 14:38
  • look at my edit, if you want to get a general solution for `x1` given `x2` the easiest way is to rearrange the ellipse matrix-form in a quadratic equation form. and that way you're getting the two `x1` solutions for each `x2` value – user2999345 Apr 06 '17 at 16:27
  • 1
    Thanks a lot! x1_1 and x2_2 have to be divided by 2A. – Astrid Apr 10 '17 at 14:55