0
$\begingroup$

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

$\endgroup$
4
  • 2
    $\begingroup$ 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). $\endgroup$
    – anderstood
    Feb 6 '21 at 18:07
  • $\begingroup$ omg! you are right! I am soo sloopy! $\endgroup$ Feb 6 '21 at 18:17
  • 2
    $\begingroup$ Terms lite d0r^(1/g) pretty much preclude obtaining an answer from NSolve. $\endgroup$
    – bbgodfrey
    Feb 7 '21 at 0:43
  • 1
    $\begingroup$ Right, that fractional power will be a problem for NSolve. Are nsolutions being sought that satisfy some criteria e.g. reals, positives? $\endgroup$ Feb 7 '21 at 16:14
2
$\begingroup$

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,
     Thread[{vars, vars /. initsolns[[1]]}]], 
   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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.