Having trouble solving a complicated system of equations

Question

I've been trying to solve a system of three equations using the 4th as a parameter, but it hasn't been working. Mathematica just seems to run indefinitely even if I set a given value for the parameter.

First, I've needed to solve a few complicated linear equations using matrix methods.

This is the matrix that I've needed to invert and I've done so without too much hassle.

z = 
  {{-1, p, q, r, (1 - p - q - r), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
   {(1 - p - q - r), -1, p, q, r, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
   {(1 - s)*r, (1 - s)*(1 - p - q - r), -1, (1 - s)*p, (1 - s)*q, 0, 0, 0, s*p, s*r, s*(1 - p - q - r), 0, 0, 0, 0, 0}, 
   {s*q, s*r, s*(1 - p - q - r), -1, s*p, (1 - s)*(1 - p - q - r), (1 - s)*q, (1 - s)*p, 0, 0, 0, 0, 0, 0, 0, 0}, 
   {p, q, r, (1 - p - q - r), -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
   {0, 0, 0, p, 0, -1, r, q, 0, 0, 0, 0, 0, 0, 0, 0}, 
   {0, 0, 0, r, 0, q, -1, (1 - p - q - r), 0, 0, 0, 0, 0, 0, 0, 0}, 
   {0, 0, 0, s*(1 - p - q - r), 0, s*r, s*p, -1, 0, 0, 0, 0, 0, (1 - s)*p, (1 - s)*q, (1 - s)*(1 - p - q - r)}, 
   {0, 0, (1 - p - q - r), 0, 0, 0, 0, 0, -1, q, r, 0, 0, 0, 0, 0}, 
   {0, 0, q, 0, 0, 0, 0, 0, r, -1, p, 0, 0, 0, 0, 0}, 
   {0, 0, (1 - s)*p, 0, 0, 0, 0,  0, (1 - s)*q, (1 - s)*(1 - p - q - r), -1, s*p, s*r,  s*(1 - p - q - r), 0, 0}, 
   {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, (1 - p - q - r), -1, q, r, 0, 0}, 
   {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, q, r, -1, p, 0, 0}, 
   {0, 0, 0, 0, 0, 0, 0, s*(1 - p - q - r), 0, 0, (1 - s)*p, (1 - s)*q, (1 - s)*(1 - p - q - r), -1, s*p, s*r}, 
   {0, 0, 0, 0, 0, 0, 0, r, 0, 0, 0, 0, 0, (1 - p - q - r), -1, q}, 
   {0, 0, 0, 0, 0, 0, 0, p, 0, 0, 0, 0, 0, q, r, -1}}

Then to get the solutions to these linear equations I've set:

a = 
  (Inverse[z]).
  ({{0}, {0}, {-s*q}, {0}, {0}, {0}, {0}, {0}, {-p}, {-(1 - p - q - r)}, {-(1 - s)*r}, {0}, {0}, {0}, {0}, {0}})
 b = 
   (Inverse[z]).
   ({{0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {-s*q}, {-p}, {-(1 - p - q - r)}, {-(1 - s)*r}, {0}, {0}})
c = (
  Inverse[z]).
  ({{0}, {0}, {0}, {0}, {0}, {0}, {0}, {-(1 - s)*r}, {0}, {0}, {0}, {0}, {0}, {-s*q}, {-p}, {-(1 - p - q - r)}})

Then I essentially want a fixed point of these equations with the top entry of a equal to p, the top entry of b equal to q and the top entry of c equal to r. If I could get this as a function of s, that would be ideal, but if that's not possible, even numerical solutions for some s would be useful.

Plugging in:

f=Part[a, 1],
g=Part[b, 1], 
h=Part[c, 1],
    NSolve[
      {f - p == 0, g - q == 0, h - r == 0, s == 0.5}, 
      {p, q, r}]

Just gives me a very long run time and seemingly no answer. Through other methods, I know that p = 0.355536, q = r = s - p = 0.5 - 0.355536 is one solution, but NSolve doesn't seem to work here and just seems to run for a long time.

Is there another way to find the solutions to these equations in terms of s or even numerically for given values of s?

I would be grateful if anyone could help or suggest anything I might be able to do to solve complicated systems of equations like mine.

Answer 1

The matrix contains an unknown parameter m (typo?). We set m = 0, then the code for finding p,q,r has the form:

m = 0; z[p_, q_, r_, 
  s_] := {{-1, p, q, r, (1 - p - q - r), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0}, {(1 - p - q - r), -1, p, q, r, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
   0}, {(1 - s)*r, (1 - s)*(1 - p - q - r), -1, (1 - s)*p, (1 - s)*q, 
   0, 0, 0, s*p, s*r, s*(1 - p - q - r), 0, 0, 0, 0, 0}, {s*q, s*r, 
   s*(1 - p - q - r), -1, 
   s*p, (1 - s)*(1 - p - q - r), (1 - s)*q, (1 - s)*p, 0, 0, 0, 0, 0, 
   0, 0, 0}, {p, q, r, (1 - p - q - r), -1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0}, {0, 0, 0, p, 0, -1, r, q, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 
   0, r, 0, q, -1, (1 - p - q - r), 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0,
    s*(1 - p - q - r), 0, s*r, s*p, -1, 0, 0, 0, 0, 
   0, (1 - s)*p, (1 - s)*q, (1 - s)*(1 - p - q - r)}, {0, 
   0, (1 - p - q - r), 0, 0, 0, 0, 0, -1, q, r, 0, 0, 0, 0, 0}, {0, 0,
    q, 0, 0, 0, 0, 0, r, -1, p, 0, 0, 0, 0, 0}, {0, 0, (1 - s)*p, 0, 
   0, 0, 0, 0, (1 - s)*q, (1 - s)*(1 - p - q - r), -1, s*p, s*r, 
   s*(1 - p - q - r), 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 
   0, (1 - p - q - r), -1, q, m, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    q, r, -1, p, 0, 0}, {0, 0, 0, 0, 0, 0, 0, s*(1 - p - q - r), 0, 
   0, (1 - s)*p, (1 - s)*q, (1 - s)*(1 - p - q - r), -1, s*p, 
   s*r}, {0, 0, 0, 0, 0, 0, 0, r, 0, 0, 0, 0, 0, (1 - p - q - r), -1, 
   q}, {0, 0, 0, 0, 0, 0, 0, p, 0, 0, 0, 0, 0, q, r, -1}}

z1[p_, q_, r_, s_] := Inverse[z[p, q, r, s]];

a[p_, q_, r_, s_] := 
  z1[p, q, r, 
    s].({{0}, {0}, {-s*
       q}, {0}, {0}, {0}, {0}, {0}, {-p}, {-(1 - p - q - 
         r)}, {-(1 - s)*r}, {0}, {0}, {0}, {0}, {0}});
b[p_, q_, r_, s_] := 
  z1[p, q, r, 
    s].({{0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {-s*
       q}, {-p}, {-(1 - p - q - r)}, {-(1 - s)*r}, {0}, {0}});
c[p_, q_, r_, s_] := 
  z1[p, q, r, 
    s].({{0}, {0}, {0}, {0}, {0}, {0}, {0}, {-(1 - s)*
       r}, {0}, {0}, {0}, {0}, {0}, {-s*q}, {-p}, {-(1 - p - q - r)}});
eq[p_, q_, r_, s_] := {Part[a[p, q, r, s], 1] - p == 0, 
   Part[b[p, q, r, s], 1] - q == 0, Part[c[p, q, r, s], 1] - r == 0};

pp[0] = 1; qq[0] = 1; rr[0] = 1; s = .5; n = 20;
Do[pp[i + 1] = First[Part[a[pp[i], qq[i], rr[i], s], 1]]; 
  qq[i + 1] = First[Part[b[pp[i], qq[i], rr[i], s], 1]]; 
  rr[i + 1] = First[Part[c[pp[i], qq[i], rr[i], s], 1]];, {i, 0, 
   n}] // AbsoluteTiming

Iterations converge quickly (0.01 sec on my comp), as can be seen from Fig. 1. The result {pp[n + 1], qq[n + 1], rr[n + 1]}={0.345055, 0.138447, 0.139816} is different than expected at s=.5, apparently m is not 0.

The cycle Do can be replaced by NestList[] or FixedPoint[] as follows

NestList[{First[Part[a[#[[1]], #[[2]], #[[3]], .5], 1]], 
   First[Part[b[#[[1]], #[[2]], #[[3]], .5], 1]], 
   First[Part[c[#[[1]], #[[2]], #[[3]], .5], 1]]} &, {1., 1., 1.}, 20]

(*Out[]= {{1., 1., 1.}, {-0.175536, -0.237042, 0.177609}, {0.178135, 
  0.110173, -0.11315}, {0.200683, 0.118494, 0.144633}, {0.284405, 
  0.126614, 0.138238}, {0.314544, 0.134986, 0.142278}, {0.331058, 
  0.136705, 0.140824}, {0.338435, 0.137736, 0.140402}, {0.341974, 
  0.138112, 0.140085}, {0.343613, 0.138295, 0.139946}, {0.344382, 
  0.138376, 0.139876}, {0.344741, 0.138414, 0.139844}, {0.344909, 
  0.138432, 0.139829}, {0.344987, 0.13844, 0.139822}, {0.345023, 
  0.138444, 0.139819}, {0.34504, 0.138446, 0.139817}, {0.345048, 
  0.138447, 0.139816}, {0.345052, 0.138447, 0.139816}, {0.345054, 
  0.138447, 0.139816}, {0.345055, 0.138447, 0.139816}, {0.345055, 
  0.138447, 0.139816}}*)

And

 With[{s = .5}, 
 FixedPoint[{First[Part[a[#[[1]], #[[2]], #[[3]], s], 1]], 
    First[Part[b[#[[1]], #[[2]], #[[3]], s], 1]], 
    First[Part[c[#[[1]], #[[2]], #[[3]], s], 1]]} &, {1., 1., 1.}]]

(*Out[]= {0.345055, 0.138447, 0.139816}*)

If m=r then we have expected result at s=.5

With[{s = .5}, 
     FixedPoint[{First[Part[a[#[[1]], #[[2]], #[[3]], s], 1]], 
        First[Part[b[#[[1]], #[[2]], #[[3]], s], 1]], 
        First[Part[c[#[[1]], #[[2]], #[[3]], s], 1]]} &, {1., 1., 1.}]]

{0.355536, 0.144464, 0.144464}

Define a function

nl[x_, n_] := 
 With[{s = x, k = n}, 
  NestList[{First[Part[a[#[[1]], #[[2]], #[[3]], s], 1]], 
     First[Part[b[#[[1]], #[[2]], #[[3]], s], 1]], 
     First[Part[c[#[[1]], #[[2]], #[[3]], s], 1]]} &, {1., 1., 1.}, 
   k]]

Using the function nl, we plot p[s],q[s],r[s] as follows

{ListLinePlot[Table[{x, nl[x, 50][[51, 1]]}, {x, .05, .95, .01}], 
  AxesLabel -> {"s", "p"}], 
 ListLinePlot[Table[{x, nl[x, 50][[51, 2]]}, {x, .05, .95, .01}], 
  AxesLabel -> {"s", "q"}], 
 ListLinePlot[Table[{x, nl[x, 50][[51, 3]]}, {x, .05, .95, .01}], 
  AxesLabel -> {"s", "r"}]}