# The solutions of system of equations are not continuous

I'm solving the system of equations that is equivalent to polynomial equation and I'm trying to plot the solutions as a function of one parameter. The problem is that instead of continuous functions I have switching between two solutions at some point. Below you can see as two solutions plotted in red and orange flips in the middle. Still it is clear that there are tree distinct branches. Further I will choose the only one branch that is always positive. But to do so I need to define the function which describe only this branch.

Here is the code to solve the system of equations:

Manipulate[
m = n + 1;
solution = Solve[{r*M0 ==
P + M0*Sum[
i*Product[Evaluate[Symbol["K" <> ToString[j]]]*1, {j, i}]*
P^i, {i, n}]/
Sum[Product[Evaluate[Symbol["K" <> ToString[j]]]*1, {j, i}]*
P^i, {i, 0, n}],
M == M0/Sum[
Product[Evaluate[Symbol["K" <> ToString[j]]]*1, {j, i}]*
P^i, {i, 0, n}],
Table[Evaluate[Symbol["MP" <> ToString[j]]]*1 ==
Times @@ (Table[
Evaluate[Symbol["K" <> ToString[i]]]*1, {i, j}]) M P^j, {j,
n}]} //
Flatten, {P, M,
Table[Evaluate[Symbol["MP" <> ToString[j]]]*1, {j, n}]} //
Flatten];, {n, 3, 1, 1}]


And here is the code to plot the solutions as a functions of parameter 'r':

valP = P /. Take[solution, {1, 4}];

molP[rVal_, K1Val_, K2Val_, K3Val_, M0Val_, DifInd_: 0] :=
Re@With[{r = rVal, K1 = K1Val, K2 = K2Val, K3 = K3Val, M0 = M0Val, índex = DifInd}, (Evaluate@D[valP, {r, índex}])];

Manipulate[
Plot[Evaluate[{molP[r, 10^k1, 10^k2, 10^k3, 10^p0, 0]}],
{r, 0.1, rrmax}, ImageSize -> 500, Frame -> True],
{{rrmax, 10, "r"}, 1, 10},
{{k1, 6, "k1"}, 2, 9},
{{k2, 3, "k2"}, 2, 9},
{{k3, 6, "k3"}, 2, 9},
{{p0, -3.5, "10^P0"}, -1, -10},
ControlPlacement -> Left]


Parameter n in Manipulate define the system of equations and at n=1 or n=2 there is no problems with branches switching:

The problems starts only when n=3. So, the question is - How I can mix the solutions in such way, so I can define all branches as continuous functions ?

Other example is below

solution = Solve[{
K1*P*L == A,
K2*P*L*L == B,
K2*A*L*L == F,
P0 == P + A + B + F,
r*P0 == L + A + 2*B + 3*F
}, {P, L, A, B, F}];

valF = F /. Take[solution, {1, 4}];

FF[rVal_, K1Val_, K2Val_, P0Val_, DifInd_: 0] :=
Re@With[{r = rVal, K1 = K1Val, K2 = K2Val, P0 = P0Val,
índex = DifInd}, (Evaluate@D[valF, {r, índex}])];

Plot[Evaluate[FF[r, 10^6, 10^7, 10^-3.8, 0]],
{r, 0.05, 2},
PlotStyle -> {Red, Blue, Green, Magenta},
PlotRange -> {{0, 1.5}, {-0.00005, 0.0001}},
ImageSize -> Large
]


I need to define 4 continuous functions instead of 4 current not continuous functions!

• Try using Solve with Quartics->False. Commented May 16, 2017 at 22:49
• Can you please fix your code so that it indeed gives an output? Just call Quit and try the Manipulate block yourself. You that there are several undefined variables including P, M0, and so on. Commented May 16, 2017 at 22:50
• @CarlWoll So easy! Thank you! Commented May 16, 2017 at 22:51
• @CarlWoll It is not universal solution ;( if I plug another parameters into molP function the solutions changes the branches even with Quartics->False option. molP[r, 10^5.34, 1000, 10^6.1, 10^-3.5, 0] Commented May 17, 2017 at 19:34
• Interesting problem. If you set all of the parameters = 1 so r is the remaining variable, then you can see internally there is a rational polynomial inside a square root where the numerator and denominator both go to zero simultaneously. At that point, the left limit is negative and the right limit is positive. Taking the absolute value gives you a smooth curve that is positive, minus the discontinuity. When it goes negative, the square root is complex. this gives the break in the curve. Commented May 19, 2017 at 22:15

For the sake of plotting, NDSolve is more likely to track a evolution of a root continuously, since it uses the derivative to predict its next value. The large number of parameters in the OP's problem lead to a correspondingly large system here. The parameters could be incorporated into a ParametricNDSolve code, but for sake of developing and debugging a working solution, it seemed convenient to have access to the components in global variables.

The following shows how to do it for the OP's "other example." Note it uses the solution from Solve to set up the initial conditions. The system algsys consists of equations for all four roots. These equations are identical except for the names of the variables. They are differentiated before being passed to NDSolve. The root each subsystem represents is determined by the initial conditions algics.

varsP = {P1, P2, P3, P4};     (* one variable for each solution *)
varsL = {L1, L2, L3, L4};
varsA = {A1, A2, A3, A4};
varsB = {B1, B2, B3, B4};
varsF = {F1, F2, F3, F4};
allvars = {{varsP, varsL, varsA, varsB, varsF}, {P, L, A, B, F}};
algsys = {K1*P*L == A, K2*P*L*L == B, K2*A*L*L == F,
P0 == P + A + B + F, r*P0 == L + A + 2*B + 3*F} /.
Equal -> Subtract /. MapThread[#2 -> Through[#1[r]] &, allvars] // Flatten;
algics = MapThread[Through[#1[r]] == (#2 /. solution) &, allvars] /. r -> 0.05;
Clear[FF2];
FF2[K1Val_, K2Val_, P0Val_, DifInd_: 0] :=
Block[{K1 = K1Val, K2 = K2Val, P0 = P0Val, índex = DifInd},
Through[varsF[r]], {r, 0.05, 2}]
];


The real and imaginary parts (green and magenta have the same real parts, blue and red have the same imaginary parts):

Plot[Evaluate[Re@FF2[10.^6, 10.^7, 10^-3.8, 0],
{r, 0.05, 2}, PlotStyle -> {Red, Blue, Green, Magenta}]


Plot[Evaluate[Im@FF2[10.^6, 10.^7, 10^-3.8, 0]],
{r, 0.05, 2}, PlotStyle -> {Red, Blue, Green, Magenta}]


• Is it possible to use FF2 as a model for NonlinearModelFit ? Commented Jul 10, 2017 at 8:37
• One could use Indexes[] to fit a particular Part[] perhaps. One should be able to fit the vector of roots Commented Jul 10, 2017 at 12:45
• @ФилиппЦветков Sorry, on a phone. Hit Send accidentally -- I won't be able to check all day, so that's really just a guess at this point. (And it's Indexed[], not Indexes[].) Commented Jul 10, 2017 at 12:49
• For instance it's continuous function of r, but what about the other parameters? Need to check Commented Jul 10, 2017 at 12:52
• Sorry! Me too.. on vocations. Cannot check (( Commented Jul 12, 2017 at 7:49

This answer uses my result from another thread, Is there a way to sort NSolve solution (roots) automatically? It doesn't use any information about the equations, it just looks to identify smooth curves from the data.

solution =
Solve[{K1*P*L == A, K2*P*L*L == B, K2*A*L*L == F,
P0 == P + A + B + F, r*P0 == L + A + 2*B + 3*F}, {P, L, A, B, F}];

valF = F /. Take[solution, {1, 4}];

FF[rVal_, K1Val_, K2Val_, P0Val_, DifInd_: 0] :=
Re@With[{r = rVal, K1 = K1Val, K2 = K2Val, P0 = P0Val,
índex = DifInd}, (Evaluate@D[valF, {r, índex}])];


Instead of immediately plotting, first generate a list of points, then manipulate them, and plot that. Here's the data.

data = Table[
Evaluate[FF[r, 10^6, 10^7, 10^-3.8, 0]], {r, 0.05, 2, .01}];

ListPlot[data // Transpose, Joined -> True]


Below is the kernel of a FoldList call. You may need to fiddle with the WEIGHT constant to get the results to come out just right. It controls how much the derivative information is used in discriminating the curves.

WEIGHT = 100;
permMatchWeight[{vals_, dels_}, l2_] :=
Module[{l1 = {vals, dels}, ps, bestPerm},
ps = Map[({#, WEIGHT (# - vals)}) &, Permutations[l2]];
bestPerm = Sort[ps, Norm[l1 - #1] < Norm[l1 - #2] &] // First
];


Here's where the data gets manipulated.

res = FoldList[permMatchWeight, {{0, 0, 0, 0}, {0, 0, 0, 0}}, data] // Rest;
orderedCurves = Map[(# // First) &, res] // Transpose;


The orderedCurvesis now a list of 4 sets of points, each tied to a single curve. On examination, you find the first two curves are identical.

ListLinePlot[orderedCurves , PlotStyle -> {Red, Blue, Green, Magenta}, ImageSize -> Large ]


In generating the data, I ignored tracking the independent variable, so you'd need to massage the result to get it back in. You could do an interpolating fit at this point to get a function for each of the curves.