# Numerically finding a curve that satisfies an underdetermined set of equations

I am trying to determine a phase diagram along the lines of Fig 3. or Fig 4. of this pubilcation. The curves plotted are binodals and are determined by the conditions of equality of chemical potentials of each component in all phases and equality of osmotic pressure in each phase (Section 4 of the SI is a good primer).

As an example, I use free energy function similar to the one defined in the paper (although the one I am actually interested in a bit more complex. Ideally I can come up with a solution that is robust and generally applicable). This is still a good starting point because my system does reduce to this under certain limits.

We narrow our attention to the case of two solutes (1 solvent) and 2 phases. Then we wish to find the locus of points in ($$\phi_a,\phi_b$$) for which the following conditions hold

$$\mu_a^{(1)} (\phi^{(1)}_a,\phi^{(1)}_b) = \mu_a^{(2)} (\phi^{(2)}_a,\phi^{(2)}_b)$$ $$\mu_b^{(1)} (\phi^{(1)}_a,\phi^{(1)}_b) = \mu_b^{(2)} (\phi^{(2)}_a,\phi^{(2)}_b)$$ $$\Pi^{(1)} (\phi^{(1)}_a,\phi^{(1)}_b) = \Pi^{(2)} (\phi^{(2)}_a,\phi^{(2)}_b)$$

Essentially, we have 4 unknowns $$\phi^{(1)}_a,\phi^{(1)}_b, \phi^{(0)}_a,\phi^{(0)}_b$$ and 3 equations. Subscripts here in my notation label the components and superscripts the phase. The solutions should look like the curves present in Fig 3 or Fig 4 of the linked publication (depending on parameter values). The conditions described above should be equivalent to equations (39)-(41) in this SI

Here are the function definitions

     F[a_, b_, uaa_, ubb_, uab_, na_, nb_] :=
a/na*Log[a] + b/nb*Log[b] + (1 - a - b)*Log[1 - a - b + 10^-14] -
1/2*uaa*a^2 - 1/2*ubb*b^2 - uab*b*a

det[a_, b_, uaa_, ubb_, uab_, na_, nb_] :=
Det[D[F[a, b, uaa, ubb, uab, na, nb], {{a, b}, 2}]] // Evaluate
\[Mu]a[a_, b_, uaa_, ubb_, uab_, na_, nb_] :=
D[F[a, b, uaa, ubb, uab, na, nb], a] // Evaluate
\[Mu]b[a_, b_, uaa_, ubb_, uab_, na_, nb_] :=
D[F[a, b, uaa, ubb, uab, na, nb], b] // Evaluate
\[CapitalPi]F[a_, b_, uaa_, ubb_, uab_, na_,
nb_] := \[Mu]a[a, b, uaa, ubb, uab, na, nb]*
a + \[Mu]b[a, b, uaa, ubb, uab, na, nb]*b -
F[a, b, uaa, ubb, uab, na, nb] // Evaluate


My first thought was to use NSolve

    NSolve[{\[Mu]a[a1, b1, 0, 0, 4.36, 10, 6] == \[Mu]a[a2, b2, 0, 0,
4.36, 10, 6],
\[Mu]b[a1, b1, 0, 0, 4.36, 10, 6] == \[Mu]b[a2, b2, 0, 0, 4.36, 10,
6],
\[CapitalPi]F[a1, b1, 0, 0, 4.36, 10, 6] == \[CapitalPi]F[a2, b2, 0,
0, 4.36, 10, 6], a1 >= 0, b1 >= 0, a2 >= 0, b2 >= 0,
a2 + b2 <= 1, a1 + b1 <= 1
}, {a1, b1, a2, b2}, Reals]


but it seems like that is not the right way to do it.

Alternatively, I thought of using FindRoot combined with grid search in the region, but I am not quite sure how to implement it, or if it is the optimum most efficient way to solve this problem. An intelligent way of choosing seed values for FindRoot might be to use the spinodal which can be visualized as follows. Since if I understand correctly the binodal curve we are interested in should envelop this so we can try seeding points on the exterior of this?

Manipulate[ContourPlot[det[a, b, uaa, ubb, uab, 10, 6] == 0,
{a, 0, 1}, {b, 0, 1}],
{uaa, 0, 10}, {uab, 0, 10}, {ubb, 0, 10}]


Q1 : Is there any other built-in function in mathematica that can help me solve my problem? Q2 : If not are there any standard approaches/algorithms that can be applied to this situation ?

EDIT : Partially Cross-posted Stack Overflow. Includes a janky-python brute force solution attempt using a grid search I attempted.

• Why not to use equations (1), (2), (3) from the paper cited? Aug 25, 2023 at 6:43
• Thank you for engaging in discussion. yes, putting $i=1,2$ in those partial derivatives and eliminating the $\eta, \pi$ terms reduces to (39)-(41) from the paper's SI (linked above) or the definitions I have here?
– jcp
Aug 25, 2023 at 10:34
• I notice there is a slight mismatch in the notation. in equations (1)-(3) in the paper, the subscript $i$ denotes the phase.
– jcp
Aug 25, 2023 at 10:37
• We can directly minimize F as well. Aug 26, 2023 at 2:52
• I mean that equations (1), (2), (3) are constraints for optimization problem for F. As I understand from discussion about data shown in Figure 2, they simulate some biological noise. Aug 26, 2023 at 16:05

We can compute data shown in Figure 4 from the paper cited using NMinimize and expression for the free energy as follows

lst = Join[{{0, 0}},
Table[{x, x Log[x]}, {x, 10^-3, 1, 10^-3}]]; xlog =
Interpolation[lst];

F[a_, b_] :=
Module[{uaa = 0, ubb = 0, uab = 4.36, na = 10, nb = 6, F, mua, mub,
PiF, det},
1/na*xlog[a] + 1/nb*xlog[b] + xlog[1 - a - b] - 1/2*uaa*a^2 -
1/2*ubb*b^2 - uab*b*a];

freeE[x_, a1_, b1_, a2_, b2_, p1_,
p2_] := {x F[a1, b1] + (1 - x) F[a2, b2], {x >= 0, x <= 1, a1 >= 0,
a2 >= 0, b1 >= 0, b2 >= 0, a1 <= 1, b1 <= 1, a2 <= 1, b2 <= 1,
x a1 + (1 - x) a2 == p1, x b1 + (1 - x) b2 == p2}};

r = Range[0, .5, .0025]; Do[p2 = .05 i;
ls[i] = Table[{a1, b1, a2, b2, x} /.
NMinimize[
freeE[x, a1, b1, a2, b2, p1, p2], {a1, b1, a2, b2, x}][[2]] //
Quiet, {p1, r}];, {i, 7}];


Visualization

Show[ListPlot[
Table[Transpose[{ls[i][[All, 1]], ls[i][[All, 2]]}], {i, 7}],
PlotRange -> All, AxesLabel -> {"\[Phi]", "\[Psi]"}],
ListPlot[
Table[Transpose[{ls[i][[All, 3]], ls[i][[All, 4]]}], {i, 7}],
PlotRange -> All]]


Same picture we compute at uab = 3.8 and combine two pictures as

Update 1. To compute data shown in Figure 3 we use code

lst = Join[{{0, 0}},
Table[{x, x Log[x]}, {x, 10^-3, 1, 10^-3}]]; xlog =
Interpolation[lst];

F[a_, b_] :=
Module[{uaa = 1.8, ubb = 0, uab = 0, na = 10, nb = 6, F, mua, mub,
PiF, det},
1/na*xlog[a] + 1/nb*xlog[b] + xlog[1 - a - b] - 1/2*uaa*a^2 -
1/2*ubb*b^2 - uab*b*a];

freeE[x_, a1_, b1_, a2_, b2_, p1_,
p2_] := {x F[a1, b1] + (1 - x) F[a2, b2], {x >= 0, x <= 1, a1 >= 0,
a2 >= 0, b1 >= 0, b2 >= 0, a1 <= 1, b1 <= 1, a2 <= 1, b2 <= 1,
x a1 + (1 - x) a2 == p1, x b1 + (1 - x) b2 == p2}};

r = Range[0, 1, .0125]; Do[p2 = If[.1 i - p1 > 0, .1 i - p1, 0];
ls[i] = Table[{a1, b1, a2, b2, x} /.
NMinimize[
freeE[x, a1, b1, a2, b2, p1, p2], {a1, b1, a2, b2, x}][[2]] //
Quiet, {p1, r}];, {i, 0, 10, 1}];


Visualization

Show[ListPlot[
Table[Transpose[{ls[i][[All, 1]], ls[i][[All, 2]]}], {i, 0, 10, 1}],
PlotRange -> All, AxesLabel -> {"\[Phi]", "\[Psi]"}],
ListPlot[
Table[Transpose[{ls[i][[All, 3]], ls[i][[All, 4]]}], {i, 0, 10, 1}],
PlotRange -> All]]


And the same picture we compute for uaa = 1.3, finally we have

Update 2. We can use equations for $$\mu_a, \mu_b, \Pi$$ as a filter for solutions computed above as follows

lst = Join[{{0, 0}},
Table[{x, x Log[x]}, {x, 10^-3, 1, 10^-3}]]; xlog =
Interpolation[lst];

fun = Module[{uaa = 1.3, ubb = 0, uab = 0, na = 10, nb = 6, F, mua,
mub, PiF, det},
F = a/na*Log[a] + b/nb*Log[b] + (1 - a - b)*Log[1 - a - b] -
1/2*uaa*a^2 - 1/2*ubb*b^2 - uab*b*a;
det = Det[D[F, {{a, b}, 2}]];
mua = D[F, a];
mub = D[F, b];
PiF = mua*a + mub*b - F; {mua, mub, PiF, det, F}];

F[a_, b_] :=
Module[{uaa = 1.3, ubb = 0, uab = 0, na = 10, nb = 6, F, mua, mub,
PiF, det},
1/na*xlog[a] + 1/nb*xlog[b] + xlog[1 - a - b] - 1/2*uaa*a^2 -
1/2*ubb*b^2 - uab*b*a];

freeE[x_, a1_, b1_, a2_, b2_, p1_,
p2_] := {x F[a1, b1] + (1 - x) F[a2, b2], {x >= 0, x <= 1, a1 >= 0,
a2 >= 0, b1 >= 0, b2 >= 0, a1 <= 1, b1 <= 1, a2 <= 1, b2 <= 1,
x a1 + (1 - x) a2 == p1, x b1 + (1 - x) b2 == p2}};

r = Range[0, 1, 1/40]; r1 = Range[0, 9/10, 1/30]; Do[
p2 = If[i - p1 > 0, i - p1, 0];
sol[i, p1] =
NMinimize[freeE[x, a1, b1, a2, b2, p1, p2], {a1, b1, a2, b2, x}][[
2]] // Quiet;, {p1, r}, {i, r1}];

eqmua1 = fun[[1]] /. {a -> a1, b -> b1}; eqmua2 =
fun[[1]] /. {a -> a2, b -> b2};

eqmua = eqmua1 - eqmua2 // Simplify;

eqmub1 = fun[[2]] /. {a -> a1, b -> b1}; eqmub2 =
fun[[2]] /. {a -> a2, b -> b2};

eqmub = eqmub1 - eqmub2 // Simplify;
eqpi1 = fun[[3]] /. {a -> a1, b -> b1}; eqpi2 =
fun[[3]] /. {a -> a2, b -> b2};

eqpi = eqpi1 - eqpi2 // Simplify

Do[ls[i] =
Table[If[(Norm[{eqmua, eqmub, eqpi}] /. sol[i, p1] /.
Indeterminate -> 1) <= 10^-3, {a1, b1, a2, b2, x} /.
sol[i, p1], Nothing], {p1, r}];, {i, r1}];


Visualization

Show[ListPlot[
Table[Transpose[{ls[i][[All, 1]], ls[i][[All, 2]]}], {i, r1}],
PlotRange -> All, AxesLabel -> {"\[Phi]", "\[Psi]"}],
ListPlot[
Table[Transpose[{ls[i][[All, 3]], ls[i][[All, 4]]}], {i, r1}],
PlotRange -> All]]


• I am grateful for your efforts. Let me make sure if I understand correctly, freeE[x_, a1_, b1_, a2_, b2_, p1_, p2_] this is the objective function for the minimisation routing, where you use this line x F[a1, b1] + (1 - x) F[a2, b2] to enforce $V_1f_1 = V_2 f_2$ ? I don't quite follow the the what the variables p1 p2 are doing? Are they the equal composition in each phase constraints?
– jcp
Aug 29, 2023 at 17:50
• Then you are taking some slices at p2 and for each slices choosing a corresponding constraint for p1 from r and then find the optimum a1,a2,b1,b2,x?
– jcp
Aug 29, 2023 at 17:55
• Any ideas on how to deal with the noisy points in your visualization and the horizontal lines that don't lie on the binodal curve? also, if this method is robust and can work for $\chi_{\phi,\phi} = 1.3, 1.8$ plots as well? Thank you for your tremendous efforts.
– jcp
Aug 29, 2023 at 17:57
• @jcp Please, note, that $x=V_1/V, 1-x=V_2/V, p_1=\bar{\phi}, p_2=\bar {\psi}$ in equations (43), (44), (45). Aug 30, 2023 at 1:31
• @jcp I don't understand your remarks about noise. Do you prefer sketch they show in Figure 4? It depends on parametrization. I used simple model with isolines p2=const. They used something like p2=const +p1. We don't know how they filtered noise. Aug 30, 2023 at 3:04