Solve for simultaneous equations

I am trying to solve the following problem, but Mathematic is timing out. Any ideas on what to do?

Subscript[z, 1] = (b^2/2 - a + Sqrt[(a - b^2/2)^2 + 2*b^2*r])/b^2

Subscript[z, 2] = (b^2/2 - a - Sqrt[(a - b^2/2)^2 + 2*b^2*r])/b^2

Solve[{x*Subscript[z, 1]*(1/(1 + n))^(Subscript[z, 1] - 1) +
y*Subscript[z, 2]*(1/(1 + n))^(Subscript[z, 2] - 1) +
l* (-((2 j)/(-a + r)) + 2/((-2 a - b^2 + r) (1 + n))) == 0,x*Subscript[z,
1] (Subscript[z, 1] - 1)*(1/(1 + n))^(Subscript[z, 1] - 2) + y*Subscript[z,
2]*(Subscript[z, 2] - 1)*(1/(1 + n))^(Subscript[z, 2] - 2) + l* (2/(-2 a - b^2 + r) ) == 0, x*Subscript[z, 1]*(m/(1 + m))^(Subscript[z, 1] - 1) +
y*Subscript[z, 2]*(m/(1 + m))^(Subscript[z, 2] - 1) +
l *(-((2 j)/(-a + r)) + (2 m)/((-2 a - b^2 + r) (1 + m))) == 0, x*Subscript[z,
1] (Subscript[z, 1] - 1)*(m/(1 + m))^(Subscript[z, 1] - 2) +
y*Subscript[z,
2]*(Subscript[z, 2] - 1)*(m/(1 + m))^(Subscript[z, 2] - 2) +
l *(2/(-2 a - b^2 + r) ) == 0}, {x, y, n, m}]


If it helps, I can assume r>0 and 0<=n,m<=1

• Are there any known constraints on the variables or parameters? Are you looking for all solutions or just real solutions? Oct 26 '20 at 22:53
• @BobHanlon I added assumptions on r and n,m to the question in case it helps. Real solutions would be ideal, but I'll take any solution. Oct 27 '20 at 6:46
• @Bill I didn't get a result from the Reduce code you shared Oct 27 '20 at 6:49

This system can be solved by noticing the first two equations do not contain $$m$$ and are linear in $$x$$ and $$y$$. So, there may be solutions for $$(x,y)=f(n)$$. Likewise, the last two equations do not contain $$n$$ and are linear in $$x$$ and $$y$$, so there may be solutions for $$(x,y)=g(m)$$. The solutions for $$(x,y)$$ will be the same, provided $$f(n)$$ and $$g(m)$$ are the same constant (independent of $$n$$ and $$m$$). That's the plan.

Rewrite the equations without subscripts

orig = {x*z1*(1/(1 + n))^(z1 - 1) + y*z2*(1/(1 + n))^(z2 - 1) +
l*(-((2 j)/(-a + r)) + 2/((-2 a - b^2 + r) (1 + n))) == 0,
x*z1 (z1 - 1)*(1/(1 + n))^(z1 - 2) +
y*z2*(z2 - 1)*(1/(1 + n))^(z2 - 2) + l*(2/(-2 a - b^2 + r)) == 0,
x*z1*(m/(1 + m))^(z1 - 1) + y*z2*(m/(1 + m))^(z2 - 1) +
l*(-((2 j)/(-a + r)) + (2 m)/((-2 a - b^2 + r) (1 + m))) == 0,
x*z1 (z1 - 1)*(m/(1 + m))^(z1 - 2) +
y*z2*(z2 - 1)*(m/(1 + m))^(z2 - 2) + l*(2/(-2 a - b^2 + r)) ==
0};


where

zrules = {z1 -> (b^2/2 - a + Sqrt[(a - b^2/2)^2 + 2*b^2*r])/b^2,
z2 -> (b^2/2 - a - Sqrt[(a - b^2/2)^2 + 2*b^2*r])/b^2}


Introduce constants

To make the equations a little more wieldy, introduce constants $$c, c_0, c_n, c_m$$.

crules = {c -> l (-((2 j)/(-a + r))),
cn -> (2 l)/ (-2 a - b^2 + r),
c0 -> (2 l)/(-2 a - b^2 + r),
cm -> (2 l)/(-2 a - b^2 + r)};

eqns = {
c + cn/(1 + n)  +  (1/(1 + n))^(-1 + z1) x z1  +
(1/(1 + n))^(-1 + z2) y z2   ==   0,
c0  +              (1/(1 + n))^(-2 + z1) x (-1 + z1) z1 +
(1/(1 + n))^(-2 + z2) y (-1 + z2) z2   ==   0,
c + cm*m/(1 + m) + (m/(1 + m))^(-1 + z1) x z1  +
(m/(1 + m))^(-1 + z2) y z2   ==   0,
c0  +              (m/(1 + m))^(-2 + z1) x (-1 + z1) z1  +
(m/(1 + m))^(-2 + z2) y (-1 + z2) z2  ==  0}


Check

Check that the original equations can be recovered by applying crules to eqns:

orig[[All, 1]] == eqns[[All, 1]] /. crules // Simplify
(*  True  *)


Solve the linear equations

Solve the first 2 equations for $$(x,y)$$ in terms of variable $$n$$

sol1 = Solve[eqns[[1 ;; 2]], {x, y}] // Flatten

(*  {x -> ((1/(1 + n))^(2 - z1) (-c - c0 - cn - c n + c z2 + cn z2 +
c n z2))/(z1 (z1 - z2)),
y -> ((1/(1 + n))^(2 - z2) (-c - c0 - cn - c n + c z1 + cn z1 +
c n z1))/(z2 (-z1 + z2))}  *)


Solve the last 2 equations for $$(x,y)$$ in terms of variable $$m$$

sol2 = Solve[eqns[[3 ;; 4]], {x, y}] // Flatten

(* {x -> (m (m/(1 + m))^-z1 (-c - c m - c0 m - cm m + c z2 + c m z2 +
cm m z2))/((1 + m)^2 z1 (z1 - z2)),
y -> (m (m/(1 + m))^-z2 (-c - c m - c0 m - cm m + c z1 + c m z1 +
cm m z1))/((1 + m)^2 z2 (-z1 + z2))  *)


Solve for $$n$$ and $$m$$

The first solution $$(x,y)$$ contains an expression for $$x$$ in terms of $$n$$. We know from the second solution that $$x$$ is independent of $$n$$, so the expression in the first solution must be a constant. Likewise for the expression for $$y$$. Here we set either derivative, $$x'(n)$$ or $$y'(n)$$, from the first solution equal to zero, solve for $$n$$, and apply our replacement rules to get a solution for $$n$$

soln = Solve[D[x /. sol1, n] == 0, n] /. crules /. zrules // Simplify // Flatten;
soln = Solve[D[y /. sol1, n] == 0, n] /. crules /. zrules // Simplify // Flatten

(*  {n -> -1 + 1/j}  *)


Now use the second solution for $$(x,y)$$, take the derivatives wrt $$m$$ and solve for $$m$$. It's a little tricky because we end up with 3 distinct solutions for $$m$$. The solutions to the equations $$x'(m)=0$$ and $$y'(m)=0$$ can be combined an simplified like this:

Solve[D[x /. sol2, m] == 0, m] /. crules // Flatten;
Solve[D[y /. sol2, m] == 0, m] /. crules // Flatten;

solm = Union[%, %%] // Sort;
solm[[1]] = First[solm] /. zrules // Simplify;

solm    (*  {m -> j/(1-j),
m -> 0        if Re[z1] < -1,
m -> 0        if Re[z2] < -1}  *)


Solve for $$x$$ and $$y$$

Now we can use the solution for $$n$$, soln, to obtain $$(x,y)$$ from sol1, or we can use one of the solutions for $$m$$, solm[[1]], to obtain the same expressions for $$(x,y)$$ from sol2.

solxy = sol2 /. solm[[1]] /. crules // Simplify

(*  {x -> (2 j^(2 - z1) l (-r + b^2 (-1 + z2) + a z2))  /
((a - r) (2 a + b^2 - r) z1 (z1 - z2)) ,
y -> (2 j^(2 - z2) l (-r + b^2 (-1 + z1) + a z1))  /
((a - r) (2 a + b^2 - r) z2 (-z1 + z2)) }  *)


Verification

We can verify the solution by applying our 3 replacement rules to the original equations, like this

orig /. solxy /. soln /. First[solm] // Simplify

(*  True, True, True, True  *)