Solving this nonlinear system of equations

Question

Suppose I have this set of equations:

$$a = x + z\qquad (1)$$ $$b = y + \frac{z}{2}\qquad (2)$$ $$ z = k_0x\sqrt{y}\qquad (3)$$

Where $a$, $b$ $\in \mathbb{R}$ and $k_0 > 0$. The values of $a$ and $b$ assume that neither x nor y can be less than zero, how do I robustly solve this system? I have attempted the Newton-Raphson method and got the following residual and Jacobian:

$$R(x,y) = \begin{Bmatrix}x+z-a \\y+\frac{z}{2}-b \\ k_0x\sqrt{y} - z\end{Bmatrix} $$ $$J(x,y) = \begin{bmatrix}1 && 0 && 1 \\ 0 && 1 && \frac{1}{2} \\ \frac{k_0\sqrt{y}}{2} && \frac{k_0x}{2\sqrt{y}} && -1 \end{bmatrix}$$

and here's the MATLAB code:

% User-defined inputs
a = 1;
b = 0;
k0 = 10;

% Initial guess
%u = [1;1;1];
u = rand(3,1)*2;

% Newton-Raphson solver
i = 0; max_iters = 30; TOL = 1e-12;
fprintf('a = %f\nb = %f\nk0 = %f\nInitial guess: (%f,%f,%f)\n\n\tResidual_x\tResidual_y\tResidual_z\tResidual_rel\n',a,b,k0,u(1),u(2),u(3));
while (true)

    % Residual function
    F = [u(1)+u(3)-a;u(2)+0.5*u(3)-b;k0*u(1)*sqrt(u(2))-u(3)];

    % Initial residual
    if i == 0
        R_init = norm(F);
    % Maximum iterations
    elseif i == max_iters
        fprintf('\tSolution did not converge.\n');
        break;
    end

    % Absolute residual
    R = norm(F);

    % Relative residual
    R_rel = R/R_init;
    fprintf('\t%e\t%e\t%e\t%e)\n',norm(F(1)),norm(F(2)),norm(F(3)),R_rel);

    % Check for convergence
    if (R_rel < TOL)
        fprintf('\tConverged.\n');
        break;
    end

    % Jacobian matrix
    J = [1,0,1;0,1,0.5;u(1)*k0*sqrt(u(2)),0.5*k0*u(1)/sqrt(u(2)),-1];

    % Solve and update
    u = u-J\F;
    i = i + 1;
end
fprintf('\nF

it seems this does not always work. If for instance at iterate $i$ I have $y^i = 0$, the solver blows up (because $\frac{1}{\sqrt{0}}$ results in an error). It seems the initial guesses also seems to play a significant role because if you run the above code with randomized guesses you will get a different solution each time.

That said, how else (or how better) can I solve what I have above? Or if I stick with this Newton-Raphson scheme, how should I choose my initial guesses?

Answer 1

This can be easily reduced to a sequence of one-variable equations with the help of sage, or any other CAS. (I replaced $\sqrt{y}\mapsto y$ to simplify things, the only thing that changes is that below only solutions with $y\geq0$ would be valid.)

R = PolynomialRing(QQ, 'a,b,k,x,y,z', order='invlex')
a,b,k,x,y,z = R.gens()
I = R.ideal(x+z-a, y*y+z/2-b, k*x*y-z)
# should be 3, meaning finitely many solutions for each set of parameters
print(I.dimension()) 
I.groebner_basis()

$$ \begin{aligned} \big[&z + x - a, \\&y^{2} - \tfrac{1}{2} x - b + \tfrac{1}{2} a, \\&x y - a y + \tfrac{1}{2} k x^{2} + b k x - \frac{1}{2} a k x, \\&a k y - \frac{1}{2}k^{2} x^{2} - b k^{2} x + \tfrac{1}{2} a k^{2} x + x - a, \\&k^{2} x^{3} + 2 b k^{2} x^{2} - a k^{2} x^{2} - 2 x^{2} + 4 a x - 2 a^{2}\big] \end{aligned} $$ Thus $x$ has to satisfy the cubic equation $$ -2 a^2+4 a x-(2+a k^2-2 b k^2) x^2+k^2 x^3 = 0 $$ (The other command to try is I.elimination_ideal([y,z]).)

Depending on $a,k$, there are at most three possible values for $x$. For each possible value of $x$, you can then solve $a=x+z$, followed by the linear equation $y=z/(kx)$, giving at most three distinct solutions. Each solution can then be checked for validity and discarded if necessary.