# 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. – KraZug Mar 9 '18 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 '18 at 14:52
• @KraZug Thanks for the comment. I have revised the post. Please see my update. – user55777 Mar 10 '18 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! – user55777 Mar 12 '18 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 – Akku14 Mar 12 '18 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. – user55777 Mar 13 '18 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. – Akku14 Mar 13 '18 at 9:14