1
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Newton's method indeed only shows local convergence (i.e., if you start close enough to the root, where "close enough" is hard to quantify). You could try a (preconditioned) fixed-point iteration, or method that are tailored to systems of multivariable polynomials (if you take the square of (3)), e.g., based on Emiris' method or homotopy methods. You might also find a symbolic solution using Gröbner bases. $\endgroup$ – Christian Clason Oct 4 '15 at 22:35
  • $\begingroup$ For given a,b,k0, sympy (using Buchberger's algorithm) spits out the full list of solutions in a second, so even if you have to solve this polyinomial system repeatedly, a symbolic approach might be worth it. $\endgroup$ – Christian Clason Oct 4 '15 at 22:49
  • $\begingroup$ The code for the J matrix and the math equation seem different for the J(3,1) element, is this no problem? Also, if you substitute (1) and (2) into (3), it appears to result in a cubic equation for z, which is solvable analytically? $\endgroup$ – roygvib Oct 5 '15 at 0:22
2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.