How to get real solutions of a "system of nonlinear equations"?

Hi Matlab community, so currently I'm writing a Matlab project to find Max and Min of an equation under a constraint using Lagrange method.

It is fine even if you don't know what Lagrange method is, because my only problem, like the title, is that: I can't find real solutions by solving system of non-linear equations.

Here's my code:

clc
clear
close all

syms x
syms y
syms lambda

%equation
z1 = 2*x^2 -2*x*y+y^3;
%constraint
z2 = x^2 + y^2 -4 == 0;
%Lagrange equation
L = z1 + lambda*lhs(z2);
%It's fine if you don't understand the code above, because the problem
%starts from here
%Here is the system of non-linear equations that I have to solve
%The three variables of this system are x, y and lambda.
Lx = diff(L,x) == 0;
Ly = diff(L,y) == 0;
Ll = diff(L,lambda) == 0;

eqns = [Lx; Ly; Ll];
vars = [x y lambda];

%Find solutions of Lx, Ly, Ll
[solx soly soll] = solve(eqns,vars);
double([solx soly soll])

Here's the output:

ans =

   1.4419 + 0.0000i   1.3860 + 0.0000i  -1.0387 + 0.0000i
   1.9092 + 0.0000i  -0.5956 + 0.0000i  -2.3120 + 0.0000i
  -2.1557 + 0.6776i   1.1704 + 1.2481i  -2.3285 - 0.6822i
  -2.1557 - 0.6776i   1.1704 - 1.2481i  -2.3285 + 0.6822i
  -1.3243 + 0.0000i   1.4987 + 0.0000i  -3.1317 + 0.0000i
  -0.3820 + 0.0000i  -1.9632 + 0.0000i   3.1394 + 0.0000i

So the problem is: I just need real solutions, which are the 1st, 2nd, 5th, 6th. So how can I get rid of the non-real solutions, which are the 3rd and 4th?

I have tried this way:

  %Find solutions of Lx, Ly, Ll
[solx soly soll] = solve(eqns,vars,'Real', true);
double([solx soly soll])

But it ended up with this error:

    Warning: Solutions are parameterized by the symbols: u. To include
    parameters and conditions in the solution, specify the 'ReturnConditions'
    value as 'true'. 
    > In sym/solve>warnIfParams (line 475)
      In sym/solve (line 364)
      In b2b (line 24) 
    Warning: Solutions are only valid under certain conditions. To include
    parameters and conditions in the solution, specify the 'ReturnConditions'
    value as 'true'. 
    > In sym/solve>warnIfParams (line 478)
      In sym/solve (line 364)
      In b2b (line 24) 
    Error using symengine
    Unable to convert expression into double array.
    
    Error in sym/double (line 698)
            Xstr = mupadmex('symobj::double', S.s, 0);
    
    Error in b2b (line 25)
    double([solx soly soll])

If you know how to fix this, please send me help, I really appreciate that.

Thanks for helping me.

Answer 1

[Alan Stevens]

If you are willing to do the differentiation by hand then the following works

% L = 2*x^2 -2*x*y+y^3 + lambda*(x^2 + y^2 -4);
% dLdx = 4*x -2*y + 2*lambda*x
% dLdy = -2*x +3*y^2 +2*lambda*y
% dLdlambda = x^2 + y^2 -4;
% Let u = [x; y; lambda];

dLdx = @(u) norm([4*u(1)-2*u(2)+2*u(3)*u(1);
             -2*u(1)+3*u(2)^2+2*u(3)*u(2);
             u(1)^2+u(2)^2-4]);

x0 = 1; y0 = -1; lambda0 = -1;
U0 = [1 1 1;
      -1 1 1;
       1 -1 -1;
      -1 -1 -1];
  [r,c] = size(U0);
  u = zeros(r,3);
  Fval = zeros(r,1);
options = optimset('Tolx',10^-12);
for i = 1:r
    u0 = U0(i,:);
    [u(i,:), Fval(i)] = fminsearch(dLdx, u0,options);
end

disp(u)

    1.4419    1.3860   -1.0387
   -0.3820   -1.9632    3.1394
    1.9092   -0.5956   -2.3120
   -1.3243    1.4987   -3.1317

disp(Fval)

   1.0e-11 *

    0.1908
    0.1591
    0.3614
    0.2479