# Having trouble solving a complicated system of equations

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.

• The code contains typos. Formulate the math problem first. – Alex Trounev Mar 7 at 12:42

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


• 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! – Mateyman Mar 9 at 12:14
• What exact result do you want to get? – Alex Trounev Mar 9 at 16:41
• 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) – Mateyman Mar 9 at 18:58
• How do you want to use this solution? – Alex Trounev Mar 9 at 23:38
• 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. – Mateyman Mar 9 at 23:48