NSolve a system of equations with a parameter

Question

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?

Answer 1

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}}} *)

Answer 2

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}
]