Is there any predictor-corrector method in Mathematica for solving nonlinear system of algebraic equations?
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!
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'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.