# Solving A System of Two Implicit Equations

I am solving a system of 6 equations, and I am interested to know the solution for different rhoBar values. Here is the input parameters:

(*Input Parameters*)
kB = 1.380649*10^-23;
absTemp = 298.15;
beta = 1/(kB*absTemp);
eps0 = 8.8541878128*10^(-12);
eps = 78.6;
f = 96485.33212;
e = 1.60217663*10^(-19);
nA = 6.02214076*10^(23);
u2 = -8/beta;
u1 = -8/beta;
u3 = -8/beta;
vW = 30*10^(-30);
vB = 5*10^(-30);
vD = 32*10^(-30);
vBnb = 5*vBb;
vBb = 30*10^(-30);
vCb = 100*10^(-30);
dLayBb = 250*10^(-12);
dLayBnb = 250*10^(-12);
dLayCb = 250*10^(-12);
k0 = 10^(-6.8)*nA;
pSite = 4.8*10^18;
cToBRatio = 0;
dToBRatio = 2;
c = -3;
rhoBar = 1.823348001*nA/dToBRatio;


The equations are:

eqn1Squared = (beta*e^2)/(
2*rhoBar*eps*eps0)*(-2*pBnb - 2*pBb + pC + pSiO)^2 ==(*(10^(-8)*
nA^2)/(dToBRatio*rhoBar^2) Exp[beta*(-e*phi0 +vD*
osmP0)]*)+Exp[beta*(-2*e*phi0 - vB*osmP0)] +
cToBRatio*Exp[beta*(e*phi0 - vC*osmP0)] +
dToBRatio*Exp[beta*(e*phi0 - vD*osmP0)] + c;
eqn2 = pSiO ==
pSite*(k0*dToBRatio*rhoBar*Exp[beta*(e*phi0 - vD*osmP0)])/(
10^(-8)*nA^2 +
k0*dToBRatio*rhoBar*Exp[beta*(e*phi0 - vD*osmP0)]);
eqn3 = (pBnb - pC)/(pSiO - pBnb - 2*pBb) ==
Exp[-beta*(u1 + 2*e*phi0 + osmP0*vBnb)];
eqn4 = (Sqrt[2]*pBb)/(pSiO - pBnb - 2*pBb) ==
Exp[-0.5*beta*(u3 + 2*e*phi0 + osmP0*vBb)];
eqn5 = pC/(pBnb - pC) == Exp[beta*(u1 - u2 + e*phi0 - osmP0*vCb)];
eqn6 = 1 - pBb*vBb/dLayBb - pBnb*vBnb/dLayBnb - pC*vCb/dLayCb -
rhoBar*Exp[beta*(-2*e*phi0 - vB*osmP0)]*vB -
cToBRatio*rhoBar*Exp[beta*(e*phi0 - vC*osmP0)]*vC -
dToBRatio*rhoBar*Exp[beta*(e*phi0 - vD*osmP0)]*
vD == (1 - rhoBar*vB - cToBRatio*rhoBar*vC -
dToBRatio*rhoBar*vD)*Exp[-beta*vW*osmP0];


Solving eqn2, eqn3, eqn4, eqn5, I managed to decouple the system of 6 equations into system of 2 equations, with variable phi0 and osmP0 left to be solved.

solution2345 = Solve[{eqn2, eqn3, eqn4, eqn5}, {pSiO, pBnb, pBb, pC}];


I then substitute solutions from eqn2 to eqn5 into eqn1Squared and eqn6

eqn1Squared = eqn1Squared /. solution2345[[1]];
eqn1Squared = eqn1Squared /. {osmP0 -> osmP0*10^6};
eqn1Squared = eqn1Squared /. {phi0 -> phi0*10^-3};
funcPhi0[phi0_, osmP0_] = (eqn1Squared[[1]] - eqn1Squared[[2]]);

eqn6 = eqn6 /. solution2345[[1]];
eqn6 = eqn6 /. {osmP0 -> osmP0*10^6};
eqn6 = eqn6 /. {phi0 -> phi0*10^-3};
funcOsmP0[phi0_, osmP0_] = (eqn6[[1]] - eqn6[[2]]);


I had tried NSolve the solve eqn6 and eqn1Squared. This did not work as it took forever to run. Then, I tried FindRoot. This method is not very reliable as I had to input guessed values for phi0 and osmP0 everytime I change the input parameter, rhoBar. To know the guessed values, I had to visually inspect the intersection between eqn6 and eqn1Squared using ContourPlot:

plot = ContourPlot[{funcPhi0[phi0, osmP0] == 0,
funcOsmP0[phi0, osmP0] == 0}, {phi0, -200, 500}, {osmP0, -2000,
2000}, MaxRecursion -> 4, PlotPoints -> 75]


I am now thinking to convert eqn6 (or its function funcOsmP0==0) into an interpolating function, osmP0 that is expressed as a function of phi0. This is so that I could substitute the function of osmP0 into eqn1Squared and solve it easily. I came across this post, and its solution could only turn the multi-line contour plot into a parametric curve.

My questions are:

1. How can I adapt its solution into my problem, and convert the implicit equation plotted in ContourPlot (f(x,y)=0) into a one to one function, or even into a one to many function (As a simple example: f(x,y)=x^2+y^2=4)?
2. Is there a more reliable and straightforward method to solve my problem? I want to solve my system of equations for various combination of input parameters, hence I am aiming to find an automated solving procedure.
• Are you looking for the intersectionpoints in your plot? Nov 11, 2023 at 21:26

First plot the function you want tointerpolate

plot2 = ContourPlot[ funcOsmP0[phi0, osmP0] == 0 , {phi0, -200, 500},{osmP0, -2000, 2000}, MaxRecursion -> 4, PlotPoints -> 75]


interpolation of plotted points

pts2=plot2[[1, 1]][[1]];
ip2 = Interpolation[pts2]
Plot[funcPhi0[phi0, ip2[phi0] ] , {phi0, ip2["Domain"][[1, 1]],
ip2["Domain"][[1, 2]]}, PlotRange -> {-1 , 1 }]


Plot shows zeros near phi0=44 and phi0=-110

NMinimize[{funcPhi0[phi0, ip2[phi0] ]^2, phi0 > 44}, phi0 ]
(*{1.93291*10^-27, {phi0 -> 44.6898}}*)

NMinimize[{funcPhi0[phi0, ip2[phi0] ]^2, phi0 < 44 }, phi0 ]
(*{6.49295*10^-14, {phi0 -> 43.6525}}*)

NMinimize[{funcPhi0[phi0, ip2[phi0] ]^2, -120 < phi0 < -100 }, phi0 ]
(*{5.6955*10^-23, {phi0 -> -110.572}}*)


Hope it helps!

Without "guessing" the phi0 try

zero1 = NMinimize[{funcPhi0[phi0, ip2[phi0]]^2 }, phi0]
(*{8.17923*10^-26, {phi0 -> 44.6898}}*)

max1 = NMaximize[{(funcPhi0[phi0, ip2[phi0]]/(phi0 - (phi0 /. zero1[[2]])))^2,ip2["Domain"][[1, 1]] < phi0 < ip2["Domain"][[1, 2]]},phi0]
zero2 = NMinimize[1/max1[[1]] ( funcPhi0[phi0, ip2[phi0]]/(phi0 - (phi0 /. zero1[[2]])))^2  , phi0]
(*{7.26201*10^-19, {phi0 -> 43.6525}}*)

max2 = NMaximize[{(funcPhi0[phi0,ip2[phi0]]/((phi0 - (phi0 /.zero1[[2]])) (phi0 - (phi0 /. zero2[[2]]))))^2,ip2["Domain"][[1, 1]] < phi0 < ip2["Domain"][[1, 2]]}, phi0]
zero3 = NMinimize[1/max2[[1]] (funcPhi0[phi0,ip2[phi0]]/((phi0 -(phi0 /. zero1[[2]])) (phi0 - (phi0 /.zero2[[2]]))))^2  , phi0]
(*{1.14371*10^-34, {phi0 -> -110.572}}*)


Stop iterated search if minimum found takes finite value !=0

• Thanks!, is there a way to solve for phi0 without requiring guessed input? Nov 12, 2023 at 8:20
• Perhaps yes, see my mofified answer Nov 12, 2023 at 9:21

An alternative approach is to use FindRoot or FindMinimum without constructing an InterpolatimgFunction, if rough estimates of the desired solutions can be provided.

FindRoot[{funcPhi0[phi0, osmP0], funcOsmP0[phi0, osmP0]}, {{phi0, -100}, {osmP0, 100}}]
(* {phi0 -> -110.572, osmP0 -> 89.7373} *)

FindMinimum[{(funcPhi0[phi0, osmP0] - funcOsmP0[phi0, osmP0])^2,
funcPhi0[phi0, osmP0] == 0, funcOsmP0[phi0, osmP0] == 0},
{{phi0, -100}, {osmP0, 100}}]
(* {2.57693*10^-15, {phi0 -> -110.572, osmP0 -> 89.7373}} *)

FindRoot[{funcPhi0[phi0, osmP0], funcOsmP0[phi0, osmP0]}, {{phi0, 50}, {osmP0, 150}}]
(* {phi0 -> 44.8119, osmP0 -> 124.638} *)

FindMinimum[{(funcPhi0[phi0, osmP0] - funcOsmP0[phi0, osmP0])^2,
funcPhi0[phi0, osmP0] == 0, funcOsmP0[phi0, osmP0] == 0},
{{phi0, 50}, {osmP0, 150}}]
(* {2.62667*10^-21, {phi0 -> 44.8119, osmP0 -> 124.638}} *)

FindRoot[{funcPhi0[phi0, osmP0], funcOsmP0[phi0, osmP0]}, {{phi0, 40}, {osmP0, 120}}]
(* {phi0 -> 43.7612, osmP0 -> 125.286} *)

FindMinimum[{(funcPhi0[phi0, osmP0] - funcOsmP0[phi0, osmP0])^2,
funcPhi0[phi0, osmP0] == 0, funcOsmP0[phi0, osmP0] == 0},
{{phi0, 40}, {osmP0, 120}}]
(* {2.91018*10^-22, {phi0 -> 43.7612, osmP0 -> 125.286}} *)


Results are identical to about ten significant figures. Which of FindRoot and FindMinimum is preferable depends on details of the problem to be solved.