# Is there any predictor-corrector method in Mathematica for solving nonlinear system of algebraic equations?

The FindRoot function in Mathematica can easily be used to solve systems of nonlinear algebraic equations. But, I want to solve a system of nonlinear equations with variations of some parameters.

However, in most cases it fails to converge. There are several methods like the arc-length method in which two parameters vary simultaneously. (See e.g. this paper)

I cannot solve my equations with NSolve, and its convergence really depends on the initial guess. So, is there any built-in function to handle this situation? For example,

ans = {x, y} /.
FindRoot[{x^2 + y^2 - #, x y - 24}, {{x, 4}, {y, 9}},
MaxIterations -> 5000] & /@ Range[1, 100, 1];
ListPlot[Table[{ans[[#, i]], Range[1, 100, 1][[#]]} & /@
Range[Length[Range[1, 100, 1]]], {i, 1, 2}], Joined -> True,
PlotRange -> All]

• If you post a specific example in Mathematica format there is a better chance of getting a useful response. Jan 21, 2016 at 23:23
• @DanielLichtblau, I found a simple example, but in fact I have a large system of nonlinear equations. Jan 22, 2016 at 8:52
• I should mention that the example given is fine for NSolve. That said, I'm not sure if your statement was in reference to this particular example or to a more difficult family of problems. Jan 22, 2016 at 16:05
• @DanielLichtblau, Thank you for your answer, these family of problems are related in tracing equilibrium of structures in finite element method. So, most of the time, there are a very large system of nonlinear equations, in which Newton method fails to converge due to limit points. Jan 22, 2016 at 17:06
• Maybe some useful info here: demonstrations.wolfram.com/… Jan 22, 2016 at 20:21

Since @hesam asked about a command, and to get a better understanding of @DanielLichtblau's approach, I tried to generalize it and package it in a function. Feedback would be appreciated!

TrackRoot[eqns_List,unks_List,{par_,parmin_?NumericQ,parmax_?NumericQ},ipar_?NumericQ,
iguess_List,opts___?OptionQ]:=

Module[{findrootopts,ndsolveopts,subrule,isol,ics,deqns,sol},
(* options *)
findrootopts=Evaluate[FindRootOpts/.Flatten[{opts,Options[TrackRoot]}]];
ndsolveopts=Evaluate[NDSolveOpts/.Flatten[{opts,Options[TrackRoot]}]];

subrule=Table[unk->unk[par],{unk,unks}];

(* use FindRoot to improve initial guess *)
isol=FindRoot[eqns/.par->ipar,Transpose[{unks,iguess}],Evaluate[Sequence@@findrootopts]];
ics=Table[{unk[ipar]==(unk/.isol)},{unk,unks}];

(* differentiate eqns *)
deqns=Map[#==0&,D[eqns/.subrule,par]];

(* track root with NDSolve *)
sol=NDSolve[Join[deqns,ics],unks,{par,parmin,parmax},Evaluate[Sequence@@ndsolveopts]][[1]];

Return[sol];
];
Options[TrackRoot]={FindRootOpts->{},NDSolveOpts->{}};

TrackRoot::usage="TrackRoot[eqns,unks,{par,parmin,parmax},initpar,initguess]
tracks a root of eqns, varying par from parmin to parmax, with initial guess
initguess at par=initpar.";


Here's an example:

{pmin, pmax} = {48.0001, 200};
tr = TrackRoot[{x^2 + y^2 - p, x y - 24}, {x, y}, {p, pmin, pmax}, 100, {10, 2}];
Plot[Evaluate[{x[p], y[p]} /. tr], {p, pmin, pmax}]


Next step I might try is to incorporate an arc-length method that could go around corners.

EDIT: Here's an attempt at using a pseudo-arclength method, inspired by this demonstration. It uses the same syntax as above and hides the fact that it uses arclength s internally. To handle multiple roots for a given parameter value it breaks the resulting InterpolatingFunction into segments between turning points (that part is particularly ugly code). Seems to work OK but I haven't tested it extensively.

TrackRootPAL[eqns_List,unks_List,{par_,parmin_?NumericQ,parmax_?NumericQ},ipar_?NumericQ,iguess_List,opts___?OptionQ]:=
Module[{s,findrootopts,ndsolveopts,smin,smax,s1,s2,subrule,isol,ics,deqns,breaks,sol,res,respart},

(* options *)
findrootopts=Evaluate[FindRootOpts/.Flatten[{opts,Options[TrackRootPAL]}]];
ndsolveopts=Evaluate[NDSolveOpts/.Flatten[{opts,Options[TrackRootPAL]}]];
smin=Evaluate[SMin/.Flatten[{opts,Options[TrackRootPAL]}]];
smax=Evaluate[SMax/.Flatten[{opts,Options[TrackRootPAL]}]];

subrule=Append[Table[unk->unk[s],{unk,unks}],par->par[s]];

(* use FindRoot to improve initial guess *)
isol=FindRoot[eqns/.par->ipar,Transpose[{unks,iguess}],Evaluate[Sequence@@findrootopts]];
ics=Join[Table[unk[0]==(unk/.isol),{unk,unks}],{par[0]==ipar}];

(* setup eqns *)
deqns=Join[
Map[#==0&,eqns/.subrule],
{Total[D[unks/.subrule,s]]^2+D[par[s],s]^2==1},
Table[unk'[0]==0,{unk,unks}],
{par'[0]==1}
];

(* track root with NDSolve *)
breaks={}; (* capture turning points *)
sol=NDSolve[Join[deqns,ics,
{WhenEvent[par'[s]==0,AppendTo[breaks,s]],
WhenEvent[par[s]==parmax,"StopIntegration"],WhenEvent[par[s]==parmin,"StopIntegration"]}],
Append[unks,par],{s,smin,smax},Evaluate[Sequence@@ndsolveopts]][[1]];

(* extract s endpoints *)
{s1,s2}=(par/.sol)["Domain"][[1]];

(* add endpoints to breaks *)
breaks=Sort[Join[{s1,s2},breaks]];

(* construct interpolatingfunctions (unk vs par) for each segment (between breaks) *)
res={};
Do[
respart={};
Do[
pts=Transpose[{(par/.sol)["Coordinates"][[1]],(par/.sol)["ValuesOnGrid"],(unk/.sol)["ValuesOnGrid"]}];
AppendTo[respart,unk->Interpolation[Select[pts,breaks[[i]]<=#[[1]]<=breaks[[i+1]]&][[All,2;;3]],
"ExtrapolationHandler"->{Indeterminate&,"WarningMessage"->False}]];
,{unk,unks}];
AppendTo[res,respart];
,{i,Length[breaks]-1}];

Return[res];
];

Options[TrackRootPAL]={FindRootOpts->{},NDSolveOpts->{},SMin->-100,SMax->100};


An example:

λ = 2.5;
tr = TrackRootPAL[{-z^3 + λ z + μ}, {z}, {μ, -10, 10}, 0, {1.4}];
Plot[Evaluate[z[μ] /. tr], {μ, -10, 10}]


Again, suggestions and improvements would be great!

• it is a great job, and as you said, maybe it is a good idea to develop it to track solutions after Limit points, Bifurcation points and etc. However, there are several advanced programs such as XPPAUT, which can be used to track solutions of nonlinear ODE, etc. Link Jan 23, 2016 at 19:37
• Thanks. I know about XPPAUT and Content, but wouldn't it be nice to have some of that functionality in Mathematica? I think so! Jan 23, 2016 at 23:02
• Is this method appropriate for an ODE with a parameter, for which I want to track the solution by varying the parameter? How do you apply boundary conditions in the code? Could you provide an example :) Thank you in advance!
– jsxs
Jun 21, 2019 at 5:31
• @jsxs It can track the equilibria of an ODE that you get by setting time derivatives equal to zero. Jun 21, 2019 at 6:12
• Hi, @ChrisK what is ipar in the argument list of TrackRoot? Should it be any value between pmin and pmax? In my example, the parameter is the domain size L. Can TrackRoot handle it? Thank you!
– jsxs
Jul 22, 2022 at 5:32

As noted by @ChrisK, this works better starting at the top. Reason being there are no real solutions below the parameter value of 48.

Using FoldList one can readily use the prior result to seed the next, that is, providing a starting point. This is a fairly common homotopy approach.

points =
Rest[FoldList[({x, y} /.
FindRoot[{x^2 + y^2 - #2, x*y - 24}, {x, #1[[1]]}, {y, #1[[2]]},
MaxIterations -> 5000]) &, {10, 5}, Range[100, 48, -1]]]

(* Out[312]= {{9.68831380576, 2.47721125483}, {9.63289204076,
2.49146361222}, {9.57705689273, 2.50598908086}, {9.5207972894,
2.5207972894}, {9.46410161514, 2.53589838486}, {9.40695767175,
2.55130307135}, {9.34935263547, 2.56702265234}, {9.29127300977,
2.58306907727}, {9.23270457346, 2.59945499274}, {9.17363232343,
2.61619379912}, {9.11404041144, 2.63329971303}, {9.05391207408,
2.65078783664}, {8.99322955501, 2.66867423468}, {8.93197401851,
2.68697602011}, {8.87012545288, 2.70571144991}, {8.80766256248,
2.72490003219}, {8.74456264654, 2.74456264654}, {8.68080146268,
2.76472167958}, {8.61635307292, 2.78540117807}, {8.55118966907,
2.80662702253}, {8.48528137424, 2.82842712475}, {8.41859601621,
2.85083165338}, {8.35109886769, 2.87387329264}, {8.28275234732,
2.89758754018}, {8.21351567389, 2.92201305177}, {8.14334446456,
2.94719204185}, {8.07219026539, 2.9731707518}, {8.,
3.}, {7.92671531783, 3.02773583227}, {7.85227181897,
3.05644029566}, {7.77659812551, 3.08618236569}, {7.69961476067,
3.11703906572}, {7.62123278463, 3.14909682963}, {7.54135211915,
3.18245317561}, {7.45985946958, 3.21721878246}, {7.37662571918,
3.25352009356}, {7.29150262213, 3.29150262213}, {7.20431854953,
3.33133520332}, {7.11487293424, 3.37321554746}, {7.02292889219,
3.41737761672}, {6.92820323028, 3.46410161514}, {6.83035261157,
3.51372782122}, {6.72895390058, 3.56667624041}, {6.62347538298,
3.62347538298}, {6.51323307597, 3.68480595122}, {6.39732143808,
3.75157012701}, {6.27449734057, 3.82500759779}, {6.14297179931,
3.90690382181}, {6., 4.}, {5.84096258932,
4.10891178175}, {5.65685424949, 4.24264068712}, {5.4244289009,
4.4244289009}, {4.89897954515, 4.89897942598}} *)


Call the parameter p. Suppose you have solved at p=100, and you want a solution at p=50. You might instead set up as a set of differential equations in a new parameter t, letting the solution set morph from the starting system at p=100 to the final one at p=50. The idea is to treat x and y as functions of t and differentiate, using the known solution at p=100 as initial values. I will skip the derivation and just show the code for this particular case.

{xsol, ysol} =
NDSolveValue[{2*x[t]*x'[t] + 2*y[t]*y'[t] == 50,
x[t]*y'[t] + y[t]*x'[t] == 0, x[0] == 9.68831380576221,
y[0] == 2.4772112548342307}, {x, y}, {t, 1, 0}]


Note that the ODE above can be handled even if we alter to let the final parameter value go below 47. To address this one would want to change the formulation into a differential algebraic system so as to enforce the algebraic equations. One way to do this is to use Method->Projection in NDSolve. I will defer to the documentation for details.

For a more complicated example, I show some path tracking code in the appendix of this work

• I'm digging this. One question: should your initial value be x[1] and y[1] rather than x[0] and y[0] (since t=0 seems to correspond to p=50 and t=1 to p=100)? I've tried both and the initial condition at t=1 version seems to match the FindRoot approach. Jan 22, 2016 at 17:52
• For me it's a bit easier to understand by sticking with p rather than t. That is {xsol, ysol} = NDSolveValue[{2*x[p]*x'[p] + 2*y[p]*y'[p] == 1, x[p]*y'[p] + y[p]*x'[p] == 0, x[100] == 9.68831380576221, y[100] == 2.4772112548342307}, {x, y}, {p, 50, 100}]; Jan 22, 2016 at 17:55
• @ChrisK Yes, I got the initial points wrong. As for using a separate parameter t, I guess it was not really needed here. Jan 22, 2016 at 19:50
• @ChrisK One reason to use a new parameter (not in this example though) is when two or more are being varied. A new parameter makes it straightforward to employ a linear homotopy between start values and end values of all the others at once. Jan 23, 2016 at 20:28
• @Daniel Lichtblau, can your method trace the solution of an ODE as a parameter varies?
– jsxs
Jul 26, 2022 at 9:07

I'd love to see a good answer to this, because it's a common problem I face. My crude improvement on your technique is to use the previous parameter value's answer as an initial guess, which helps FindRoot. Even better, use linear or quadratic extrapolation and add some adaptive step size. But you won't be going around any bends like the vertical limit point in your diagram!

ans = {};
{x0, y0} = {2.477, 9.688};
Do[
{x0, y0} = {x, y} /.
FindRoot[{x^2 + y^2 - p, x y - 24}, {{x, x0}, {y, y0}}];
AppendTo[ans, {x0, y0}];
, {p, 100, 1, -1}];
ListPlot[Table[{ans[[#, i]], Range[1, 100, 1][[#]]} & /@
Range[Length[Range[1, 100, 1]]], {i, 1, 2}], Joined -> True,
PlotRange -> All]


Notice I'm going from p=100 down to p=1, which runs more smoothly.

• What problems do you solve that require this sort of analysis? I solve cold plasma dispersion relations with numerical effects included in order to improve numerical stability of particle-in-cell codes. The typical computation involves numerous roots, some complex, that move about, often merging in frequency space as parameters are varied. I know roughly where in phase space they lie, but computations must be efficient, because parameter space is large. Jan 24, 2016 at 18:37
• @bbgodfrey For me, equilibria of nonlinear differential equation models in theoretical ecology & evolutionary biology. Jan 24, 2016 at 19:55
• @ChrisK You should check my answer for an explicit pseudo-arclength continuation method. With some effort, it can be extended to all kinds of operators. Mar 3, 2017 at 19:06

Here is my approach using an explicit pseudo-arclength method.

The typical predictor-corrector method uses the parameter as continuation parameter. For example, let $$G(u,h) = 0. \label{eq:root}\tag{1}$$ If we know a solution $$(u_0,h_0)$$, then we predict the solution at $$h_1 = h_0 + \Delta h$$ by noting that, if the Jacobian matrix $$G_u$$ is invertible, then $$u'(h_0) = - G_u^{-1}(u(h_0),h_0)G_h(u(h_0),h_0) = u_0'.$$ Then we propose the predictor $$u_p = u_0 + u_0' \Delta h.$$ and use some method (usually Newton's) to obtain the corrector.

The problem with this method is that it fails in a fold. To circumvent this problem, instead of parameterizing the solution curve by the continuation parameter, we do it by its arclength. In this case, you have to think that everything is a function of $$s$$, i.e., $$X(s) = \big(u(s),h(s)\big)$$, like so (taken from [1]):

In this scheme, not only $$u$$ is unknown, but also $$h$$. There are several ways to close the system, the most common being to use the scalar normalization $$(X_1 - X_0)^T \dot{X_0} = \Delta s.$$ This is the equation of a plane, which is perpendicular to the tangent $$\dot{X}(s)$$ at a distance $$\Delta s$$ from $$X_0$$. This plane will intersect the solution curve if $$\Delta s$$ and the curvature of the curve is not to large. So, we extend \eqref{eq:root} to $$G(u,h) = 0,$$ $$(u - u_0)^T \dot{u}_0 + (h - h_0)\dot{h}_0 - \Delta s = 0,$$ where $$(u_0,h_0) = \big(u(s_0),h(s_0)\big)$$ is a known solution. Now we only need to calculate $$\dot{u}_0$$ and $$\dot{h}_0$$ to get things going. This can be done by differentiating $$G$$ with respect to $$s$$ and using the normalization condition $$\dot{u}_0^T \dot{u}_0 + \dot{h}_0^2 = 1$$:

\begin{align}\label{eq:extroot}\tag{2} \dot{u}_0 &= -G_u^{-1}(u_0,h_0) G_h(u_0,h_0) \dot{h}_0,\\ \dot{h}_0 &= \pm\left(1+\|G_u^{-1}(u_0,h_0) G_h(u_0,h_0)\|^2\right)^{-1/2}, \end{align}

where we choose the sign depending on the direction we want to go. In this scheme, given a solution $$(u_0,h_0)$$, our predictor is \begin{align} u_p &= u_0 + \dot{u}_0 \Delta s,\\ h_p &= h_0 + \dot{h}_0 \Delta s. \end{align}

Finally, a practical way to approximate the derivatives of $$u$$ and $$h$$ is, instead of solving \eqref{eq:extroot} at each step and determining the sign, approximate this equation by $$\pmatrix{G_u(u_1,h_1) & G_h(u_1,h_1) \\ \dot{u}_0 & \dot{h}_0}\pmatrix{\dot{u}_1 \\ \dot{h}_1} = \pmatrix{0 \\ 1},$$ which has the advantage of choosing the right direction at each step.

Details on why this scheme works on folds can be found in Keller's classic notes Lectures on Numerical Methods In Bifurcation Problems.

## Some Mathematica examples:

I'm not an MMA expert; these codes are for illustrative purposes only, and haven't been crafted for speed or elegance. I'd love to see some of our house names take a shot into packing TrackRoot[..., Method -> PseudoArcLength] ;)

### OP's own example:

funG = {x^2 + y^2 - h, x y - 24};
extfunG = {x^2 + y^2 - h, x y - 24, (x - xPrev) dxPrev + (y - yPrev) dyPrev +
(h - hPrev) dhPrev - ds};
chi = LinearSolve[D[funG, {{x, y}}], -D[funG, h]];
extchi = LinearSolve[Join[D[funG, {{x, y, h}}], {{dxPrev, dyPrev, dhPrev}}], {0, 0, 1}];

{x0, y0, h0} = {x, y, 100} /. FindRoot[funG /. h -> 100, {x, 4}, {y, 9}, MaxIterations->500];

dh0 = -(1 + chi.chi)^(-1/2) /. {x -> x0, y -> y0, h -> h0}; (* left continuation *)
{dx0, dy0} = -chi dh0 /. {x -> x0, y -> y0, h -> h0};

ClearAll[predcorr];
predcorr[xi_, dxi_, step_] := Module[{xp, yp, hp},
{xp, yp, hp} = xi + step dxi;
{{x, y, h}, extchi /. {dxPrev->dxi[[1]], dyPrev->dxi[[2]], dhPrev->dxi[[3]]}, step} /.
FindRoot[extfunG /. {xPrev -> xi[[1]], yPrev -> xi[[2]], hPrev -> xi[[3]],
dxPrev -> dxi[[1]], dyPrev -> dxi[[2]], dhPrev -> dxi[[3]], ds -> step},
{x, xp}, {y, yp}, {h, hp}, MaxIterations -> 500]
]


We track the root using a NestWhileList. Note the Check function, ready to diminish the step-size in order to ensure convergence (There might be a more elegant way to address this problem).

solCurve = NestWhileList[
Check[predcorr[#1, #2, #3], {#1, #2, #3/2}] & @@ # &,
{{x0, y0, h0}, {dx0, dy0, dh0}, 0.1}, #[[1, 3]] < 101 &
] // DeleteDuplicates[#, (#1[[1]] == #1[[2]] &)] &;


Here is a bifurcation diagram:

ListPlot[{#1[[3]], #1[[1]]} & @@@ solCurve, AxesLabel -> {h, x}]


### Chris K example:

funG1 = -z^3 + 5/2 z + h;
extfunG1 = {-z^3 + 5/2 z + h, (z - zPrev) dzPrev + (h - hPrev) dhPrev - ds};
chi1 = -D[funG1, h]/D[funG1, z];
extchi1 = LinearSolve[Join[{D[funG1, {{z, h}}]}, {{dzPrev, dhPrev}}], {0, 1}];

{z0, h0} = {z, -10} /. FindRoot[funG1 /. h -> -10, {z, -3}, MaxIterations -> 500];
dh0 = (1 + chi1^2)^(-1/2) /. {z -> z0, h -> h0}; (* right continuation *)
dz0 = -chi1 dh0 /. {z -> z0, h -> h0};

ClearAll[predcorr1];
predcorr1[zi_, dzi_, step_] := Module[{zp, hp},
{zp, hp} = zi + step dzi;
{{z, h}, extchi1 /. {dzPrev -> dzi[[1]], dhPrev -> dzi[[2]]}, step} /.
FindRoot[extfunG1 /. {zPrev -> zi[[1]], hPrev -> zi[[2]], dzPrev -> dzi[[1]],
dhPrev -> dzi[[2]], ds -> step}, {z, zp}, {h, hp}, MaxIterations -> 500]
]


We track the root:

solCurve1 = NestWhileList[
Check[predcorr1[#1, #2, #3], {#1, #2, #3/2}] & @@ # &, {{z0, h0},
{dz0, dh0}, 0.1}, #[[1, 2]] < 10 &] // DeleteDuplicates[#, (#1[[1]] == #1[[2]] &)] &;


Here is the bifurcation diagram:

ListPlot[Reverse /@ solCurve1[[All, 1]], AxesLabel -> {h, z}]


• I should point out that my answer also includes a pseudo-arclength method. Jul 6, 2019 at 10:37
• @Pragabhava, can this method trace the solution of an ODE as a parameter varies?
– jsxs
Jul 26, 2022 at 8:56
• @jsxs Yes, it can; you only need to define the derivatives of $G$ correctly (i.e. using Frechet's definition). For more details, see Keller's notes on the subject. In practice, as far as MMA goes (and I recall), you can wrap your ode into a black box function using NDSolve and use the method directly as stated. See my comment here for details. Jul 26, 2022 at 23:01

I used the Mathematica built-in function NDSolve to apply the arclength continuation technique. This is what it gives for OP's example:

eqns = {x[s]^2 + y[s]^2 - p[s] == 0, x[s] y[s] - 24 == 0};
unkows = {x[s], y[s], p[s]};
norml = {D[unkows, s].D[unkows, s] - 1 == 0};
IC = {x[0] == 2, y[0] == 2, p[0] == 110, p'[0] == -1};
SolCont = NDSolve[Join[eqns, norml, IC], unkows, {s, 0, 100}];
ParametricPlot[{p[s],x[s]} /. SolCont, {s, 0, 100}]


I tested the code for very large systems, it worked well too.