3
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ The code contains typos. Formulate the math problem first. $\endgroup$ Mar 7, 2020 at 12:42

1 Answer 1

5
$\begingroup$

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. Figure 1

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

Figure 1

$\endgroup$
8
  • $\begingroup$ Ah thanks very much! m should be r. Have managed to find numerical solutions now, still can't solve exactly but thank you very much! $\endgroup$
    – Mateyman
    Mar 9, 2020 at 12:14
  • $\begingroup$ What exact result do you want to get? $\endgroup$ Mar 9, 2020 at 16:41
  • $\begingroup$ I've been trying to express all of p, q and r in terms of s (on the interval 0<=p,q,r<=1 (with 0<=p+q+r<=1 as well) $\endgroup$
    – Mateyman
    Mar 9, 2020 at 18:58
  • $\begingroup$ How do you want to use this solution? $\endgroup$ Mar 9, 2020 at 23:38
  • $\begingroup$ It would just be quite nice to get the explicit form of it so I could make surther calculations based on varying the parameter s and see how this affects p,g and r. I'd also quite like to plot p as a function of s on the interval [0,1] which I haven't been able to do currently. $\endgroup$
    – Mateyman
    Mar 9, 2020 at 23:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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