NSolve a system of equations with a parameter

I have a system of equations with a single parameter Ee that I would like to vary from 0 to 1 with step sizes of 0.01. I have been able to get solutions with this format by manually changing Ee (there's no solution for a few values of Ee but that's fine), but now I want to get these solutions automatically. This is my code right now:

T = 298
e = 1.6*10^-19
F = 96500
k = 1.38*10^-23
R = 8.31
Γ = 10^-10
A = 0
B = F/(2.3 R T)
c = 88.5
d = 7.2
H = 2 k T/e
G = -F Γ*10^6

(* Ee = ???? *)

x4 = Log10[x/(1 - x)] == A + y B
x5 = z == c w
x6 = z + v == d H y /H
x8 = Ee == w + y
x9 = v == G x

eqns = {x4, x5, x6, x8, x9}
vars = {x, y, z, w, v};
NSolve[eqns, vars, Reals]


I saw another answer around here that suggested using this:

parameter = {1,2,3,4,5};
equation = x + parameter;

NSolve[# == 0, {x}] & /@ equation


or this:

NSolve[x + # == 0, x] & /@ {1, 2, 3, 4}


But I'm not sure how to apply that to my system. I tried to replace Ee in equation x8 with # and then writing NSolve[eqns, vars, Reals] & /@ {0, 1}, but that's incorrect. I also tried setting Ee = {0,1} and writing NSolve[#, vars, Reals] & /@ eqns, but that is incorrect as well since Mathematica returns the warning "Infinite solution set has dimension at least 1." What will allow me to solve this system at every value of the parameter?

• have you tested that you get a solution for a particular Ee? You should give values for your constants. – george2079 Sep 6 '17 at 19:14
• I've been able to go through and manually enter values for Ee across my whole range and get the solutions. Now I just want to do that automatically. I'll add the constants, I guess people can't run this without them. – ham Sep 6 '17 at 19:17

You want to get the solution as a function of Ee, so define a function that given a value of Ee gives you a solution for the rest of the parameters:

solution[evalue_] := NSolve[eqns /. Ee -> evalue, vars, Reals]


The expression eqns /. Ee -> evalue takes the equations stored in eqns and replaces Ee by whatever you feed into the function (evalue). Once you've done this, you can use this function in a table just like any other function:

Table[solution[e], {e, 0, 1, 0.25}]

(* {{{x -> 0.262543, y -> -0.0264738, z -> 2.34293, w -> 0.0264738, v -> -2.53354}},
{{x -> 0.993994, y -> 0.130961, z -> 10.535, w -> 0.119039, v -> -9.59205}},
{{x -> 0.999999, y -> 0.361547, z -> 12.2531, w -> 0.138453, v -> -9.64999}},
{{x -> 1., y -> 0.592738, z -> 13.9177, w -> 0.157262, v -> -9.65}},
{{x -> 1., y -> 0.823929, z -> 15.5823, w -> 0.176071, v -> -9.65}}} *)

• That said, I think that some of these answers may not be accurate; applying these solutions to your initial equations yields False in several cases, particularly for the first equation. You may want to ask another question about how to address this if this problem persists. – Michael Seifert Sep 6 '17 at 19:58
• I think that's due to the approximations that Mathematica makes when doing NSolve. I get this warning every time: "NSolve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result." What I've been doing is feeding the results back into expressions for the left and right sides and judging if they're "close enough" by myself. This system that I'm solving is already an approximation (x6 is supposed to be a Sinh), but I'll watch out if there's a significant difference. – ham Sep 6 '17 at 20:10

I would convert your equations into an ODE and let NDSolve construct interpolating functions of the solutions as a function of the parameter Ee. Converting your equations eqns into an ODE can be done as follows:

eq2 = eqns /. Thread[vars -> Through[vars[Ee]]];
ode = D[eq2, Ee]


{((1 - x[Ee]) (x'[Ee]/(1 - x[Ee]) + ( x[Ee] x'[Ee])/(1 - x[Ee])^2))/(Log[10] x[Ee]) == 16.9427 y'[Ee], z'[Ee] == 88.5 w'[Ee], v'[Ee] + z'[Ee] == 7.2 y'[Ee], 1 == w'[Ee] + y'[Ee], v'[Ee] == -(193/20) x'[Ee]}

We also need an initial point, which can be obtained using NSolve as you did:

ip = Equal @@@ First@NSolve[eq2 /. Ee->0, Through[vars[0]], Reals]


NSolve::ratnz: NSolve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result.

{x[0] == 0.262543, y[0] == -0.0264738, z[0] == 2.34293, w[0] == 0.0264738, v[0] == -2.53354}

Now we are ready to use NDSolveValue:

{xsol, ysol,zsol,wsol,vsol}=NDSolveValue[Flatten@{ode, ip}, vars, {Ee, 0, 2}];


Visualization:

Plot[
{xsol[t], ysol[t], zsol[t], wsol[t], vsol[t]},
{t, 0, 2}
]