Up until now I have been finding the roots for coupled nonlinear equations by first plotting the curves and finding the roots for different values of parameters to find the seed (if there is no good guess Mathematica throws a warning about inability to find a sufficient decrease in merit function) and if the curves are smooth enough then first having one guess which I know from plots or analytics, and then using that to find root and then using the solution of this as a guess for the next and so on.

Question: Is there a way to automate this process of finding the guess?

The specific problem at hand is that I am trying to solve the following system of nonlinear equations :

\[Gamma][\[Mu]_, \[Beta]_, w_] := 
 1/2 + 1/\[Pi] Im[ 
      1/2 + \[Beta]/(2 \[Pi]) (I \[Mu] + w)] - (\[Beta] w)/(2 \[Pi])
       PolyGamma[1, 1/2 + \[Beta]/(2 \[Pi]) (I \[Mu] + w)]]

dlta[d_, mu_, bta_, w_, Ex_] := 
 0.5 Ex (\[Gamma][mu + d, bta, w] - \[Gamma][mu - d, bta, w])

nn1[d_, mu_, bta_, w_, Ex_] := 
 0.5 Ex (\[Gamma][mu + d, bta, w] + \[Gamma][mu - d, bta, w])

rootCell2[bta_, w_, Ex_, mu0_] := 
 FindRoot[{Abs[Abs[dlta[d, mu1, bta, w, Ex]] - Abs[d]] == 0, 
   Abs[Abs[nn1[d, mu1, bta, w, Ex]] + mu1 - mu0] == 0}, {{mu1, 1}, {d,

The trouble I am facing is that for FindRoot to converge it must have a good seed guess but I am trying to solve the above equations for a range of parameters where bta$\in[0.01,5]$, w$\in[0,5]$,mu0 $\in[-5,5]$ and Ex$\in[0.5,5]$.

Edit: The method of continuity (suggested below by userrandrand) does indeed give a guess but because of the sharp changes in the curves as a function of the parameters, I need to have tiny changes in the parameters, which makes the code slow.

Yesterday I saw that there are some answers on this site which use advantage of contour plots of one of the equations and checking sign changes on the second function as a guess. But in cases where the function is not exactly zero but close to zero this method does not work. Perhaps I should try to adapt them for some small values and not zero. I should also say that NMinimize fails to solve the system above.

  • 1
    $\begingroup$ I think information is missing that makes it difficult to answer this question. There are solutions with d==0, but probably you are not interested in them? You did not specify the range of Ex? Perhaps you can provide a concrete set of values for the arguments where FindRoot fails? Since FindRoot uses a Newton method by default, I would remove unnecessary Abs, for example your first question has the form $||a|-|b|| = 0$ which can be written (assuming $a,b$ are real) in the form $a=b$ or $a=-b$ or $a^2 = b^2$. $\endgroup$
    – user293787
    Commented Aug 11, 2022 at 20:52
  • $\begingroup$ If you have sharp analytical changes, then the NDSolve method might work as the step size used by the solver adapts to maintain accuracy. That said you should probably try to remove the Abs if it is possible. If you prefer to stay with continuity then you can make a program that recursively divides the stepsize by two each time you receive an error from FindRoot. I might have used Check and Quiet to check if FindRoot failed, I am not sure. $\endgroup$ Commented Aug 11, 2022 at 22:50
  • $\begingroup$ I made such a program (that uses only FindRoot, continuity and dividing step size by two recursively) but I usually prefer the NDSolve method where the stepsize is chosen automatically and if I remember well it is generically quicker although in some cases taking somewhat large stepsizes with FindRoot can be faster. $\endgroup$ Commented Aug 11, 2022 at 22:52
  • $\begingroup$ That said you have a rather high dimensional space of parameters such that high accuracy methods might take too much time. In that case, possible ideas might include training a neural network (personally i only tried 1 dimensional data and digit recognition) or maybe some Monte-Carlo methods not sure. $\endgroup$ Commented Aug 11, 2022 at 22:58
  • $\begingroup$ Maybe a multi dimensional Taylor expansion around some point in parameter space would shift the problem to equations on the coefficients of the expansion and would remove some of the parameters, I am not sure. That said, it restricts the domain of validity and might require considering multiple expansion points. $\endgroup$ Commented Aug 11, 2022 at 23:25

3 Answers 3


Edit: For problems with few roots in the domain of parameters considered (or low dimensional simple manifolds) and when function evaluation is not costly, it might be better to skip the long discussion below and use the visual method by @Akku14 and the numerical method by @AlexTrounev.

In the following I will discuss general methods with simple examples and in the end I will discuss a possible solution to the original problem

The question does not seem specific to mathematica and is more of a question on numerical methods. Hence, you might get a broader range of answers by looking at other parts of stack exchange that deal with numerical methods. Not sure which one, maybe the stack exchange on mathematics as there might be applied mathematicians there.

I actually worked on a similar situation and I am planning on publishing a paper that will use the method that I explain below. The solution was inspired from references I read on the homotopy analysis method (you do not need to know it to understand the following) although the explanation you will find on Wikipedia about this method might not give you the broad idea of the method.

Let $x$ be the solution of a system of potentially non linear equations denoted as ${\cal N}$ with parameters $\theta=(\theta_1,\theta_2,...)$:

${\cal N}[x,\theta]=0$

$x$ could have many components but in the examples of this first part we will focus on a 1-dimensional problem. The original problem is 2-dimensional and I will give an answer for that problem in the second part.


(*Below ${\cal N}$=eq, g=$x$ , $\theta$={ll,r} *)

eq[g_,l_,r_]=Cos[r*g]+ll*Sin[g]; eq[g,ll,r]==0

According to the implicit function theorem $x$ is a function of $\theta$. We will take the hypothesis that $x$ is a smooth function of $\theta$.

First, for all the methods below you need a good starting point. Sometimes $\theta=0$ or $\theta=1$ will suffice and maybe sometimes you should use RandomReal. If there is a limit of the parameters where the problem becomes simple then you might want to consider obtaining some analytical results using perturbation theory in that limit to get a good starting point. It is possible that the limit itself plus maybe some random noise can be enough to get a good starting point for FindRoot when the parameters are close to that limit.


Solve[eq[g,1,1],g] gives all solutions one being -Pi/4. ll=1 and r=1 will be our starting point.

The methods:

By continuity (robust and precise but slow)

Consider a starting point $\theta=\theta^{(0)}$ and a grid of (multi-dimensional) points $\theta^{(i)}$ with $\theta^{(i)}-\theta^{(i+1)}$ small. Then, you can a make a program that operates as a for loop that finds the solution for $\theta^{(i)}$ with FindRoot and that uses that solution as an initial condition for the point $\theta^{(i+1)}$.


(* intialization *)


(* iteratively obtaining roots *)




(* Verifying the roots *)

MapThread[eq[#1,ll,#2] &,{groot/@Range[iterations+1],rlist}]


(* {1.1102210^-16, 0., 0., -5.5511210^-17, 0., 5.5511210^-17, -1.1213310^-14, -4.3298710^-15, -1.8318710^-15,
-7.2164510^-16, -2.7755610^-16} *)

By differentiability (quick)

In this method you turn the problem to a differential equation by taking a derivative of the equation and then using NDSolve. I like this method because the adaptive methods of NDSolve finds the step size for you and changes it automatically if it thinks it needs a smaller step size.Also in principle this should allow bigger step sizes as the the guess for the next step uses both the value and the slope of the previous step whereas the continuity method only looks at the previous value. However, in my experience, I have struggled to get bigger step sizes then what NDSolve does automatically. I also suggest checking sometimes if the results you find with this method do solve the non linear equation you have. This might also be true with the method above.

As mentioned before the solution $x$ is actually a function of the parameters $\theta$ of the equation. Hence the equation can be written more explicitly as :

${\cal N}[x(\theta),\theta]=0$

Taking the derivative of the equation above with respect to one of the parameters in $\theta$ which we will choose to be $\theta_1$, we obtain a differential equation:

$\partial_{\theta_1} x(\theta)\partial_x{\cal N}[x(\theta),\theta]+\partial_{\theta_1}{\cal N}[x(\theta),\theta]=0$



Because of the coefficient $\partial_x{\cal N}[x(\theta),\theta]$ that multiplies $\partial_{\theta_1} x(\theta)$, in my experience, the above differential equation is often stiff and if you try to solve it with NDSolve you will receive an error (this does not happen with the example provided). However, if I remember correctly Mathematica also suggests an option that allows it to work with stiff equations. I do not remember the name of the option but if you try this and do not see the suggestion by Mathematica then please write a comment below and I will search my files for the option. Once that option is included, you can solve the differential equation with a starting point $(\theta_1^{(\textrm{start})},\theta_2^{(0)},...)$, an end point $(\theta_1^{(\textrm{end})},\theta_2^{(0)},...)$ and an initial condition for $x$ at $\theta_1^{(\textrm{start})}$ which will give you $x(\theta)$ in that range.



Testing the accuracy:


output: verification

The precision is good but not as precise as with FindRoot, this is maybe related to the automatic values in PrecisionGoal and AccuracyGoal in NDSolve. You can try to choose your own precision if you want.

In your problem this would be one line in the cube of parameters you have and you can consider the cube of points as a cube of lines in the different parameters. You can decide how you want to space the lines.

The initial condition for $\theta_1=\theta_1^{(\textrm{start})}$ can be obtained using FindRoot or analytical results in some limit of the parameters. If I remember correctly this method works only if you have a good initial condition.

If the equation has a singularity in the domain you are searching (not necessarily the function itself it could be a derivative of the function) then the adaptive methods of NDSolve will choose smaller and smaller steps until the integration is not possible. A singularity maybe means that the program is at the boundary of the manifold defined by the set of roots of the equations.

I like to see the steps that NDSolve is taking in real time when using this differential equation method. If you are interested I can search my notebook for how I did it but it probably involves StepMonitor in NDSolve and using the function Dynamic to plot the points in real time (for my problem plotting in real time did not add much extra time but for your problem it might as it is possibly simpler).

Errors can accumulate over large domains when using the NDSolve method, so you might want to check if the solution obtained by NDSolve is really a root of the initial non linear equation. Also, you might want to polish the results (increase accuracy) by using FindRoot with the result from NDSolve as an initial condition.

Edit (! These are suggestions that I have not tried really): One could also consider maybe by piecewise analyticity using a multi-dimensional Taylor-expansion in the parameters restricted perhaps to sub-domains of the original problem. Another possibility might be by simplicity using regularized neural networks that are trained on some examples. This method is by simplicity because regularized neural networks seem to have a bias towards simple functions and so if the shape of the manifold defined by the zeroes of the equations is simple, then a neural network might be a good option in some cases (maybe high dimensional ones I am not sure).

The original problem

Edit: I had modified the $\gamma$ function in an experiment and I might have not put the factors back correctly or something because after copy pasting your function again I got new results.

The main issues in applyting the method using differentiability is the Abs and Im in the function.

I removed all the Abs and you can check a posteriori whether the signs should be different.

Concerning the Im, taking a derivative you will get an expression with Im'[...].Mathematica's derivative seems to handle cases with Im badly so I defined a derivative for this problem (I took the minimal amount of properties needed for this problem, this derivative might not generalize well to other problems):

deri[a_,r_] /; FreeQ[a,r]:= 0
deri[a_*b_,r_] /; FreeQ[a,r]:= a*deri[b,r]

I define also rules for Im so that it implements basic properties:

imrule={Im[a_ + b_] -> Im[a] + Im[b], Im[bta^m_.*y_] -> bta^m*Im[y], Im[w*y_] -> w*Im[y]}

setting up the equations from the original problem (I removed the Abs !!) :

eq1[d_, mu1_, bta_, w_, Ex_]=dlta[d, mu1, bta, w, Ex]-d//Simplify; 
eq2[d_, mu1_, bta_, w_, Ex_,mu0_]=nn1[d, mu1, bta, w, Ex]+ mu1 - mu0 //Simplify;
eqlist[d_, mu1_, bta_, w_, Ex_,mu0_]={eq1[d, mu1, bta, w, Ex],eq2[d, mu1, bta, w, Ex, mu0]};

(* derivative in the direction of w *)

dweqlist[bta_, w_, Ex_, mu0_]=deri[#, w] & /@ eqlist[d[w], mu1[w], bta, w, Ex, mu0] /.imrule /.deri->D //Simplify;

(* intialization *)


parameter values: {btai, wi, Exi, mu0i}={2.65194, 1.39173, 4.91738, 4.48104}

{mu1i,di}={mu1,d}/.FindRoot[eqlist[d, mu1, btai, wi, Exi, mu0i],{{mu1,1},{d,3}}]

output: {1.09501, 1.33762}

Notice that d is non zero although it might seem that there only roots with d=0 (* Please check this as I might have missed a term while copy-pasting or something*)

To check this you can use the graphical method by @Akku14:


You can have a different view by plotting instead equation_1^2+equation_2^2 :


The roots are encircled by the red contours

(* obtaining d and mu1 for w belonging to the interval [wi, 5] *)

(Edit: I think the error I got about complex values was because of an experiment I did and I forgot to re-initialize the d[w] and mu1[w] with FindRoot)

{dsol,mu1sol}=NDSolveValue[Thread[dweqlist[btai, w, Exi, mu0i]==0]~Join~    {mu1[wi]==mu1i,d[wi]==di},

The solver shows that it is unable to integrate pass w=2.118... and for a good reason as we will see.

(* Checking that the code works *)

Plot[eqlist[dsol[w], mu1sol[w], btai, w, Exi, mu0i],{w,wi,2.11}]  (*{equation_1,equation_2}*)

output: accuracy

Let's see why the NDSolver had trouble pass w=w_c$\simeq 2.12$ :

Plot[{dsol[w], mu1sol[w]}, {w, wi, 2.11}, PlotLabels -> Automatic]

output: d_and_mu1_vs_w

The function d goes to 0 in a singular manner that looks like sqrt if I had to guess. I show a contour plot near w_c below (I lowered the contour values to have a more detailed view):


For w> w_c there is one root and for w<w_c there are three for the above chosen parameters. This is probably a pitchfork bifurcation similar to the mean-field equation of the Ising model (easy axis magnetism).

If there are physics behind these equations w=w_c might be the onset of new physics in the model (If the discussion above is not due to the fact that I omitted the Abs function)

PS: If there are physics behind those equations I am curious as to what they are as I am not used to seeing a Gamma function in the wild in physics. Otto Hölder's theorem says that the Gamma function is a hypertranscendental function (or a transcendentally transcendental function) and so it satisfies no algebraic differential equation whose coefficients are rational functions (and algebraic initial conditions). As most physics problems seem to be linked to differential equations in this category, it felt like I would never see this function in a physics problem (outside of string theory).

  • $\begingroup$ This answer is great & can be improved by including specific examples using Wolfram Language code blocks to describe how to achieve the answer to the OP’s question! As it is, the code you do include is formatted in TeX and this is not so useful to the end user (i.e., the OP & other users on this SE). Can you, please, include some example code blocks? $\endgroup$ Commented Aug 11, 2022 at 17:45
  • $\begingroup$ Hi @CATrevillian , thank you for your suggestion, I will take a simple non linear equation as an example that could help indeed. $\endgroup$ Commented Aug 11, 2022 at 17:54
  • $\begingroup$ Sorry for the large pictures.I do not know how to choose the image size on this website. I tried choosing the image size on mathematica which kind of worked but I do not know what to do when I take screenshots (like with manipulate) $\endgroup$ Commented Aug 12, 2022 at 23:12
  • 1
    $\begingroup$ For me these integrals arise in the context of generating functions of the Dirac genus (to calculate anomalies in a certain model) $\endgroup$
    – Charlie
    Commented Aug 13, 2022 at 0:33
  • $\begingroup$ Great answer, btw. Thanks! $\endgroup$
    – Charlie
    Commented Aug 13, 2022 at 0:35

Not a full answer, but a plädoyee for graphical methods to get an overview.

   Abs[Abs[dlta[d, mu1, bta, w, Ex]] - Abs[d]], {mu1, 0, 2}, {d, -1, 
    5}, Contours -> {0, .01, .1}, ContourShading -> False, 
   ContourStyle -> Red, PlotPoints -> 50, FrameLabel -> {mu1, d}], 
   Abs[Abs[nn1[d, mu1, bta, w, Ex]] + mu1 - mu0], {mu1, 0, 2}, {d, -1,
     5}, Contours -> {0, .01, .1}, ContourShading -> False, 
   PlotPoints -> 30], ImageSize -> 300], {{bta, 1}, 1/100, 5, 
  Appearance -> "Labeled"}, {{w, 2}, 0, 5, 
  Appearance -> "Labeled"}, {{Ex, 1}, 1/2, 5, 
  Appearance -> "Labeled"}, {{mu0, 2}, -5, 5, 
  Appearance -> "Labeled"}, ContinuousAction -> False]

enter image description here

  • $\begingroup$ Nice now I understand what the OP wanted, I do not know if it is simpler to look at the contour plot of equation_1^2+equation_2^2 or two contour plots of equation_1 and equation_2 as is done here. $\endgroup$ Commented Aug 12, 2022 at 8:22
  • $\begingroup$ @userrandrand , good idea. Saves time. $\endgroup$
    – Akku14
    Commented Aug 12, 2022 at 10:35

This problem can be solved for some range of parameters with Jacobian-free Newton–Krylov method in combination with FindRoot. This is only one example compare to rootCell2

\[Gamma][\[Mu]_, \[Beta]_, w_] := 
  1/2 + 1/\[Pi] Im[
       1/2 + \[Beta]/(2 \[Pi]) (I \[Mu] + 
           w)] - (\[Beta] w)/(2 \[Pi]) PolyGamma[1, 
        1/2 + \[Beta]/(2 \[Pi]) (I \[Mu] + w)]];

dlta[d_, mu_, bta_, w_, Ex_] := 
  0.5 Ex (\[Gamma][mu + d, bta, w] - \[Gamma][mu - d, bta, w]);

nn1[d_, mu_, bta_, w_, Ex_] := 
  0.5 Ex (\[Gamma][mu + d, bta, w] + \[Gamma][mu - d, bta, w]);

rootCell2[bta_, w_, Ex_, mu0_] := 
  FindRoot[{Abs[Abs[dlta[d, mu1, bta, w, Ex]] - Abs[d]] == 0, 
    Abs[Abs[nn1[d, mu1, bta, w, Ex]] + mu1 - mu0] == 0}, {{mu1, 
     1}, {d, 3}}];

Solver definition

jvapprox[v_, F_, x_, epsilon_] := 
  Module[{dim, sum, R, per, xper, y}, dim = Length[x]; 
   sum = Sqrt[epsilon] Sum[1 + x[[i]], {i, dim}]; 
   per = If[Norm[v] > epsilon, sum/(dim Norm[v]), sum/dim ]; R = F[x];
    xper = x + per v; y = (F[xper] - R)/per; y];
jfnk[F_, x0_, epsilon_, maxiter_] := 
  Module[{x, R, per, v, dim, df, x1, jv, var}, R = F[x0]; x = x0; 
   dim = Length[x0]; v = Array[u, {dim}]; 
   sol = Array[s, {maxiter + 1}]; 
   sol[[1]] = Table[u[i] -> .1, {i, dim}]; 
   Do[jv = jvapprox[v, F, x, epsilon] - R; var = v /. sol[[i]]; 
    sol[[i + 1]] = FindRoot[jv, Table[{u[j], var[[j]]}, {j, dim}]]; 
    x = x - v /. sol[[i + 1]]; R = F[x];, {i, 1, maxiter}]; {x, R}];


SeedRandom[12345]; With[{bta = RandomReal[{.01, 5}], 
  w = RandomReal[{0, 5}], mu0 = RandomReal[{-5, 5}], 
  Ex = RandomReal[{.5, 5}]}, 
 F = Y |-> {Abs[
     Abs[dlta[Indexed[Y, 1], Indexed[Y, 2], bta, w, Ex]] - 
      Abs[Indexed[Y, 1]]], 
    Abs[Abs[nn1[Indexed[Y, 1], Indexed[Y, 2], bta, w, Ex]] + 
      Indexed[Y, 2] - mu0]}] ;
x0 = {1/10, 1/10}; jfnk[F, x0, 10^-12, 4] // AbsoluteTiming

(*Out[]= {0.0367553, {{-3.02998*10^-11, 1.22223}, {2.13368*10^-11, 

Therefore we have solution {d, mu1}={-3.02998*10^-11, 1.22223} with error of {2.13368*10^-11,3.41949*10^-14} consequently. This solution we can compare to

With[{bta = RandomReal[{.01, 5}], w = RandomReal[{0, 5}], 
   mu0 = RandomReal[{-5, 5}], Ex = RandomReal[{.5, 5}]}, 
  rootCell2[bta, w, Ex, mu0]] // AbsoluteTiming


During evaluation of In[14]:= FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances.

(* Out[]= {0.0164415, {mu1 -> 0.324484, d -> -2.52861*10^-9}}*)
  • $\begingroup$ As I understand the functions jvapprox and jfnk are problem independent, that is you can choose whatever non linear function F you want. That can be quite nice to have in a toolbox. Why is it that this method works when using FindRoot directly does not ? $\endgroup$ Commented Aug 12, 2022 at 8:30
  • 1
    $\begingroup$ I have used JFNK method in a case of singular Jacobian and similar problems when we try to compute roots with controlled residuals. Unfortunately in your case there is no unique solution due to Abs and Im functions. So, may be this method is not relevant. $\endgroup$ Commented Aug 12, 2022 at 8:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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