With scalar functions it is easy to construct a mixed Newton-Raphson-bisection algorithm so that the solution always stays inside the given bounds in which it is bracketed. What if I have a vector function which I want to finds its roots (each root depending on the other roots), but I know that the true roots have the same sign as my initial guess, and that it must stay that way in every step of the iteration procedure in order that the solution doesn't become irretrievable. Can I use this information to ensure the convergence of my solution?

To be more clear, I have a vector $\pmb{f}$ of functions, which are all the same except evaluated at different positions. Each of these component functions depend on a root vector $\pmb{x}$ comprised of all the roots at these positions. Symbolically, $\pmb{f}(\pmb{x})$ looks like $$f_i(\pmb{x})=\int_{y_i-a}^{y_i}g(y^{\prime})x(y^{\prime})dy^{\prime}+h(y_i)$$ where $x=x(y)$ is a function of position. Equivalently (?), $\pmb{x}$ can be considered as a vector solution where each of the different discrete values $y$ takes gives a component $x_i=x(y_i)$.

Now I know this is an integral equation system, but mine is actually much more complex and I'd like to solve it with Newton-Raphson if possible. The problem is that with NR alone, the solution tends to wander off in a region where some of the roots change sign, and this makes sure the solution gets forever lost.

Just to make sure, the way I proceed using only Newton-Raphson to find the roots of a vector function $\pmb{f}(\pmb{x})$ is by computing $$dx_i(\pmb{x})=\frac{f_i(\pmb{x})}{f_i^{\prime}(\pmb{x})}$$ where $$f_i^{\prime}(\pmb{x})=\frac{\partial f_i(x)}{\partial x}\Bigr|_{\pmb{x}}$$ where the boldface was dropped to show the differentiation is with respect to a variable, and then applying $$\pmb{x}_{n+1}=\pmb{x}_n-\pmb{dx}_n$$ until convergence.

I tried to use some kind of bisection method (which I know doesn't work for multidimensional problems) but I can't seem to be able to bracket a root. I asked that ${x_1}_i={x_i}_0/10$ and ${x_2}_i={x_i}_0\times 10$, calculated $\pmb{f}_1=\pmb{f}(\pmb{x}_1)$, for example. But it doesn't seem to me that such method can be generalized for a vector function.

In a way, I guess I could try a multivariate Newton-Raphson method, but usually in these cases, the vector $\pmb{x}$ is a vector of distinct variables. In my case, it's a single variable but which depends on position, hence is numerically an array. Doing partial derivatives would mean doing a derivative of that variable only at a specific position, and since I have to integrate over all positions (see below for the form of the $f_i(\pmb{x})$), it would be pretty awkward I think. Or maybe I am mistaken and it has nothing to do with Newton-Raphson. At least I know the method I explained worked for a relatively simple case, as it gave the same result as another root-finding scheme.

Is there a way to solve such a problem? I guess that if I try to solve the integral equation problem I am not guaranteed that the solution won't change sign in the process either?

  • $\begingroup$ It's similar, but it's not a system of equation with multiple, distinct variables. The function which I am trying to find its roots is a vector; but each component is the same equation but at a different geometric point. And each root is a quantity at a different geometric point. Each of these component functions depend of the whole vector root. I'll update my post as it's not so clear on this. But from what I understand, in this case it doesn't make sense to use a multivariable Newton method? $\endgroup$ Commented Aug 13, 2017 at 0:16
  • $\begingroup$ updated. I also fixed an error regarding my function f $\endgroup$ Commented Aug 13, 2017 at 0:33

1 Answer 1


I have been facing this kind of problems many years ago and my initial attempt w has to adapt subroutine rtsafe of Numerical Recipes to $n$ dimensions. It revealed to not be good because one of the variables could force many bisection steps in order that all variables stay within bounds.

Finally, what I decided to do is to turn the problem into the minimization of the norm using bound constraints. We have so many good optimizers !

Just give a try and, please, let me know about your experience.

  • $\begingroup$ That's exactly what I tried to do! You mean minimizing $\pmb{f}\cdot\pmb{f}$ in the "traditionnal" way? I know that in Numerical Recipes (2nd edition) they incorporate this minimization with multivariate Newton-Raphson, warning that simply doing the minimization alone could be problematic. Is there a specific(s) method I should try in such a case? But since it's not a "true" multivariate Newton-Raphson problem and the functions and roots are closely related, the system may be better behaved? $\endgroup$ Commented Aug 13, 2017 at 3:12
  • $\begingroup$ Well, this has been my experience. Optimization with bound constraints is really the simplest. If your opimizer is not able to handle bound constraints (I doubt this could be the case), you could use a "simple" optimizer changing variable $$x_i=a_i+\frac{b_i-a_i}{1+e^{-y_i}}$$ The $y_i$ are unbounded but you could face many problems if one of the solutions is close to one of its bounds. $\endgroup$ Commented Aug 13, 2017 at 3:23
  • $\begingroup$ I must say I am not familiar with optimization algorithms. I can only find constrained maximization algorithms (save simulated annealing) in Numerical Recipes, is there any book or on-line resource you can recommend me on the subject? $\endgroup$ Commented Aug 13, 2017 at 3:43
  • $\begingroup$ What is your computing environment ? Which programming language are you using ? Have you access to subroutine libraries ? If yes, which ones ? $\endgroup$ Commented Aug 13, 2017 at 3:56
  • $\begingroup$ I work on a GNU/Linux system and I program in FORTRAN 77. I don't think I have any access to subroutine libraries. $\endgroup$ Commented Aug 13, 2017 at 15:38

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .