# Solving simultaneous nonlinear algebra equations numerically

I want to solve three algebra equations with the three unknowns $J$, $N_i$, and $T_i$. I set the values for the parameters firstly, then try to use NSolve to solve the system, and finally, want to plot $J$ as a function of $h$ as it ranges from $0.01$ to $1$ in increments of $0.1$.

pup = 101000; Ng0 = 0; R = 8.3145; \[Rho] = 0.00138458;
molL = 40449.5; Tscale = 358; r = 1.64477; Le = 2.59541;
p0 = 5330; T0 = 292.15; \[Delta] = 1; Tg0 = 0.87; Tl0 = 0.87;
Q = 8.5*10^-3; \[Kappa] = 276; c = 0.5;

ps[Ti_] := p0/pup*Exp[molL/R*(1/T0 - 1/(Ti*Tscale))]

eq1 = J == (Tl0 - Ti)/h + c J*(Tg0 - Ti)/(Exp[(J*\[Delta])/(\[Rho]*\[Kappa])] - 1);
eq2 = Q*J == ps[Ti] - Ni/(r + (1 - r)*Ni);
eq3 = (1 - Ni)/(1 - Ng0) == Exp[(-J \[Delta])/(\[Rho]*\[Kappa]*Le)];

solnJ = NSolve[{eq1, eq2, eq3}, {J, Ti, Ni}, Reals][[1, 1]]


The solnJ gives

(*J == (0.87 - Ti)/h + (0.5 J (0.87 - Ti))/(-1 + E^(2.61681 J))*)


Here, h could be considered as a parameter. The equations thus need to be solved numerically with fixed h later.

To plot, I used

ParametricPlot[{h, solnJ[h]}, {h, 1/100, 1}, Frame -> True,
Axes -> False, PlotRange -> {{0, 1}, {0, 0.04}}, AspectRatio -> 0.5]


which gives me a blank figure. I guess the equations still have not been solved numerically. And I also tried

sol = Table[{h, Evaluate[J /. solnJ]}, {h, 1/100, 1, 1/10}]


Then, I got the following warning messages:

ReplaceAll::reps: {...} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

What can be done to avoid the warning messages and solve the equations numerically? Then plot J as a function of h.

Some findings:

1. The warning arises from J /. solnJ;

2. If I give a value for h, say, $0.5$, and use NSolve[{eq1, eq2, eq3}, {J, Ti, Ni}] instead without the restriction of Reals, Mathematica runs for a long time.

• You are currently saying that J is an unknown function in your equations by writing J[h, Ti, Ni], you need to just call it J. Also c and kappa are undefined. Mar 9, 2018 at 14:09
• Welcome to Mathematica.SE. I suggest that: 1) You take the introductory Tour 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise.
– gwr
Mar 9, 2018 at 14:52
• @KraZug Thanks for the comment. I have revised the post. Please see my update. Mar 10, 2018 at 5:55

Why not use FindRoot?

fr[k_?NumericQ] :=
FindRoot[{eq1, eq2, eq3} /. h -> k, {{Ti, 1}, {Ni, 2}, {J, 3}}]

Plot[{Ti /. fr[h], Ni /. fr[h], J /. fr[h]}, {h, 0, 1},
PlotStyle -> {Red, Green, Blue}]


• thanks for your help @Akku14. Your idea is first converting the equations to three parameter equations in $k$ then using Plot to evaluate each unknowns as a function of $h$. A key technique interests me is that fr[k] is a list of rules and contains the parameter solutions of Ti, Ni, and J, so why the replace, say, Ti /. fr[h] knows which rule should be transferred to Ti? Thank you! Mar 12, 2018 at 14:16
• I used /. h -> k only, because FindRoot does not recognize, that eq1,eq2,eq2 have h as parameter. Better style is to define eq1 as eq1[h_] ,... and call it in FindRoot as eq1[h]. So I only wanted to avoid more typing. Second if you call for example fr[1/2] you see wchich solution is for Ti Mar 12, 2018 at 17:38
• Hi @Akku14 Thanks for your reply. With eq1[h_] :=..., and FindRoot[{eq1[h], eq2[h], eq3[h]}, {{Ti, 0.1}, {Ni, 0.2}, {J, 0.3}}]. It gives me an error >...not a list of numbers with dimensions {3}. Could you please update your answer. Mar 13, 2018 at 2:17
• You should also think a litte bit for your own. If you call FindRoot[{eq1[h], eq2[h], eq3[h]}, {{Ti, 0.1}, {Ni, 0.2}, {J, 0.3}}] standing alone, h of course has no numeric value and FindRoot can't work. Define fr[h_?NumericQ] as I did or call FindRoot[{eq1[1/2], eq2[1/2], eq3[1/2]}, {{Ti, 0.1}, {Ni, 0.2}, {J, 0.3}}]  with h==1/2 for example. Mar 13, 2018 at 9:14