# 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. – Carl Woll May 16 '17 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. – halirutan May 16 '17 at 22:50
• @CarlWoll So easy! Thank you! – Филипп Цветков May 16 '17 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] – Филипп Цветков May 17 '17 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. – MikeY May 19 '17 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 ? – Филипп Цветков Jul 10 '17 at 8:37
• One could use Indexes[] to fit a particular Part[] perhaps. One should be able to fit the vector of roots – Michael E2 Jul 10 '17 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[].) – Michael E2 Jul 10 '17 at 12:49
• For instance it's continuous function of r, but what about the other parameters? Need to check – Michael E2 Jul 10 '17 at 12:52
• Sorry! Me too.. on vocations. Cannot check (( – Филипп Цветков Jul 12 '17 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.