# NSolve for system of 14 non-linear equations

I am really new to Mathematica. Trying to solve the following system of equations (below). Get an error from Mathematica and have no idea what is the problem. Will very much appreciate any help.

b = 0.99; g = 0.1119; psi = 9*10^(-5);
z0 = 1; z1 = 1; w0w = 0.0008; w = 0.46; c = 0.67; sr = 0.67; th = \
0.25; ar = 0.0074; ah = 0.03; sw = 0.88; w0r = 0.0008;
e = 0;
wh = 0.006;

NSolve[c1 - b *(1 + r1w) == 0,
d0r^(1/g) - (1/psi)*c1*1/(r1w - r1r) == 0,
k0h - (1/ah)*(1/r1w)*z1 - q0 == 0,
q0*k0w - (1 - w) bl -
w0w*(b*(1 - sw)*r1w)/(th - b*(1 - sw)*(rk1w - r1w)) == 0,
q0*k0w - w0w - d0w - bl ==
0, (q0 + ar*k0r)*k0r - c*bl -
w0r*((b*(1 - sr)*r1r)/(th - b*(1 - sr)*(rk1r - r1r))) ==
0, (q0 + ar*k0r)*k0r - w0r - d0r - bl == 0, k0w + k0r + k0h == 1,
c1 + w0w*((rk1w - r1w)*((q0*k0w - (1 - w) bl)/w0w) + r1w) +
w0r*((((q0 + ar*k0r)*k0r - c*bl)/w0r) + r1r) - z1 - z0*wh == 0,
c0 - z0*wh - w0w - w0r + ah*(k0h^2)/2 + ar*(k0r^2)/2 == 0,
rk1r - r1r - (rb1 + e - r1r)/c == 0,
rk1w - r1w - (rb1 + e - r1w)/w == 0, rk1r - z1/(q0 + ar*k0r) == 0,
rk1w - z1/q0 == 0, {rb1, c1, c0, r1w, r1r, rk1w, rk1r, bl, d0r, d0w,
q0, k0w, k0r, k0h}
]


• From the documentation on NSolve, the syntax is NSolve[expr,vars]. You need to wrap your equations in a List with braces. Whether NSolve will manage to find the solution(s) is another problem. If you know some appropriate initial guesses or if you are looking for a single solution, you might help Mathematica by using FindRoot (one solution vs exhaustive search). Commented Feb 6, 2021 at 18:07
• omg! you are right! I am soo sloopy! Commented Feb 6, 2021 at 18:17
• Terms lite d0r^(1/g) pretty much preclude obtaining an answer from NSolve. Commented Feb 7, 2021 at 0:43
• Right, that fractional power will be a problem for NSolve. Are nsolutions being sought that satisfy some criteria e.g. reals, positives? Commented Feb 7, 2021 at 16:14

You can obtain some solutions as follows. First, a small change turns the system into polynomials. We use NSolve to get solutions, and use each of those as initial values for FindRoot using the original equations. It turns out that we require fairly high precision in order to guarantee small residuals, so I rationalized everything.

b = 99/100;
g = 1119/10000;
psi = 9*10^(-5);
z0 = 1;
z1 = 1;
w0w = 8/10000;
w = 46/100;
c = 67/100;
sr = 67/100;
th = 1/4;
ar = 74/10000;
ah = 3/100;
sw = 88/100;
w0r = 8/10000;
e = 0;
wh = 6/1000;

vars = {rb1, c1, c0, r1w, r1r, rk1w, rk1r, bl, d0r,
d0w, q0, k0w, k0r, k0h};


Here I change a fractional exponent in the second expression into something symbolic.

exprs0 = {c1 - b*(1 + r1w),
d0r^(1/g1) - (1/psi)*c1*1/(r1w - r1r),
k0h - (1/ah)*(1/r1w)*z1 - q0,
q0*k0w - (1 - w) bl -
w0w*(b*(1 - sw)*r1w)/(th - b*(1 - sw)*(rk1w - r1w)),
q0*k0w - w0w - d0w - bl,
(q0 + ar*k0r)*k0r - c*bl -
w0r*((b*(1 - sr)*r1r)/(th - b*(1 - sr)*(rk1r - r1r))),
(q0 + ar*k0r)*k0r - w0r - d0r - bl,
k0w + k0r + k0h - 1,
c1 + w0w*((rk1w - r1w)*((q0*k0w - (1 - w) bl)/w0w) + r1w) +
w0r*((((q0 + ar*k0r)*k0r - c*bl)/w0r) + r1r) - z1 - z0*wh,
c0 - z0*wh - w0w - w0r + ah*(k0h^2)/2 + ar*(k0r^2)/2,
rk1r - r1r - (rb1 + e - r1r)/c,
rk1w - r1w - (rb1 + e - r1w)/c,
rk1r - z1/(q0 + ar*k0r),
rk1w - z1/q0};
exprs1 = Numerator[Together[exprs0 /. g1 -> 1]];


Obtain initial solutions.

Timing[initsolns =
NSolve[exprs1, vars, WorkingPrecision -> 500];]

(* Out[760]= {8.3716, Null} *)


Now use the first of these to get solutions to the original system.

N[rt1 = FindRoot[(exprs0 /. g1 -> g) == 0,
Evaluate@Apply[Sequence,
WorkingPrecision -> 500, MaxIterations -> 500]]

(* Out[774]= {rb1 -> -1.66326, c1 -> -2.05359, c0 -> -2.35161,
r1w -> -3.07433, r1r -> -2.95291, rk1w -> -0.968261,
rk1r -> -1.02807, bl -> -11.7898, d0r -> 3.89187, d0w -> 6.87657,
q0 -> -1.03278, k0w -> 4.7565, k0r -> 8.11875, k0h -> -11.8753} *)


Check residuals:

In[776]:= exprs0 /. g1 -> g /. rt1

(* Out[776]= {0.*10^-500, 0.*10^-493, 0.*10^-499, 0.*10^-497, 0.*10^-499,
0.*10^-499, 0.*10^-499, 0.*10^-499, 0.*10^-499, 0.*10^-500,
0.*10^-499, 0.*10^-499, 0.*10^-500, 0.*10^-500} *)


Other initial solutions may not fare this well. Also I notice that the result is not close to the initial point. In a case like this one might do better to set up a homotopy between solutions of the initial system and solutions of the desired system. This can be done as a system of ODEs.