# Speeding up solving a system of simultaneous equations

I am trying to solve the following system of equations. Solve takes a very long time to solve the following system of equations. Is there any way to speed it up?

l = 4;
vars = Flatten[Table[{a[i], b[i], c[i], d[i]}, {i, l}]];

eq1 ={{x (x a[1] + x b[1] + d[1] (1 - 2 z[1]) +
c[1] (-1 + 2 z[1]))}, {x (x a[1] - x b[1] + 2 c[1] z[1] +
2 d[1] z[1])}} == {{1}, {0}};

eq2={{a[4]}, {c[4]}} == {{0}, {0}};

eq3=Table[
{{E^(x h[i]) x a[i]+E^(-x h[i]) x b[i]+E^(-x h[i]) d[i] (1+x h[i]-2 z[i])+E^(x h[i]) c[i] (-1+x h[i]+2 z[i])},{E^(x h[i]) x a[i]-E^(-x h[i]) x b[i]+E^(-x h[i]) d[i] (-x h[i]+2 z[i])+E^(x h[i]) c[i] (x h[i]+2 z[i])},{-((E^(x h[i]) x a[i] (1+z[i]))/y[i])+(E^(-x h[i]) x b[i] (1+z[i]))/y[i]+(E^(x h[i]) c[i] (2-x h[i]-4 z[i]) (1+z[i]))/y[i]+(E^(-x h[i]) d[i] (2+x h[i]-4 z[i]) (1+z[i]))/y[i]},{-((E^(x h[i]) x a[i] (1+z[i]))/y[i])-(E^(-x h[i]) x b[i] (1+z[i]))/y[i]+(E^(-x h[i]) d[i] (1-x h[i]) (1+z[i]))/y[i]-(E^(x h[i]) c[i] (1+x h[i]) (1+z[i]))/y[i]}}=={{x a[1+i]+x b[1+i]+d[1+i] (1-2 z[1+i])+c[1+i] (-1+2 z[1+i])},{x a[1+i]-x b[1+i]+2 c[1+i] z[1+i]+2 d[1+i] z[1+i]},{-((x a[1+i] (1+z[1+i]))/y[1+i])+(x b[1+i] (1+z[1+i]))/y[1+i]+(c[1+i] (2-4 z[1+i]) (1+z[1+i]))/y[1+i]+(d[1+i] (2-4 z[1+i]) (1+z[1+i]))/y[1+i]},{-((x a[1+i] (1+z[1+i]))/y[1+i])-(x b[1+i] (1+z[1+i]))/y[1+i]-(c[1+i] (1+z[1+i]))/y[1+i]+(d[1+i] (1+z[1+i]))/y[1+i]}}
,{i, l - 1}];

Solve[Join[{eq1, eq2}, eq3], vars]


The notebook file containing equations in matrix form can be downloaded from the following link: Solve Problem.nb

-
I have deleted Sjoerd's answer at his request because he stated that it was incorrect. – Mr.Wizard Dec 31 '13 at 17:36

One thing you can do is replace the exponential terms E^(-x h[i]) with a variable e[i], solve the system, and replace the exponential terms at the end. Taking away the exponentials limits what Mathematica can do with the system. If it can solve this more generic system, we still win. And since there is less that can be tried, it is likely to be faster.

Update: Polynomial equations are also easier to solve. Originally, I thought rational functions might be converted to polynomial equations, but apparently not. We can multiply eq3 through by the denominators, expand, and simplify. Since I'm changing eq3, here is the complete code.

Clear[a, b, c, d, h, y, z, l];
l = 4;
vars = Flatten[Table[{a[i], b[i], c[i], d[i]}, {i, l}]];

eq1 = {{x (x a[1] + x b[1] + d[1] (1 - 2 z[1]) +
c[1] (-1 + 2 z[1]))}, {x (x a[1] - x b[1] + 2 c[1] z[1] +
2 d[1] z[1])}} == {{1}, {0}};

eq2 = {{a[4]}, {c[4]}} == {{0}, {0}};

eq3 = Table[
y[1 + i] y[i] e[i] * (* multiply by common denominator *)
{{E^(x h[i]) x a[i] + E^(-x h[i]) x b[i] +
E^(-x h[i]) d[i] (1 + x h[i] - 2 z[i]) +
E^(x h[i]) c[i] (-1 + x h[i] + 2 z[i])}, {E^(x h[i]) x a[i] -
E^(-x h[i]) x b[i] + E^(-x h[i]) d[i] (-x h[i] + 2 z[i]) +
E^(x h[i]) c[
i] (x h[i] + 2 z[i])}, {-((E^(x h[i]) x a[i] (1 + z[i]))/
y[i]) + (E^(-x h[i]) x b[i] (1 + z[i]))/
y[i] + (E^(x h[i]) c[i] (2 - x h[i] - 4 z[i]) (1 + z[i]))/
y[i] + (E^(-x h[i]) d[i] (2 + x h[i] - 4 z[i]) (1 + z[i]))/
y[i]}, {-((E^(x h[i]) x a[i] (1 + z[i]))/
y[i]) - (E^(-x h[i]) x b[i] (1 + z[i]))/
y[i] + (E^(-x h[i]) d[i] (1 - x h[i]) (1 + z[i]))/
y[i] - (E^(x h[i]) c[i] (1 + x h[i]) (1 + z[i]))/y[i]}} ==
y[1 + i] y[i] e[i] * (* multiply by common denominator *)
{{x a[1 + i] + x b[1 + i] + d[1 + i] (1 - 2 z[1 + i]) +
c[1 + i] (-1 + 2 z[1 + i])}, {x a[1 + i] - x b[1 + i] +
2 c[1 + i] z[1 + i] +
2 d[1 + i] z[
1 + i]}, {-((x a[1 + i] (1 + z[1 + i]))/
y[1 + i]) + (x b[1 + i] (1 + z[1 + i]))/
y[1 + i] + (c[1 + i] (2 - 4 z[1 + i]) (1 + z[1 + i]))/
y[1 + i] + (d[1 + i] (2 - 4 z[1 + i]) (1 + z[1 + i]))/
y[1 + i]}, {-((x a[1 + i] (1 + z[1 + i]))/
y[1 + i]) - (x b[1 + i] (1 + z[1 + i]))/
y[1 + i] - (c[1 + i] (1 + z[1 + i]))/
y[1 + i] + (d[1 + i] (1 + z[1 + i]))/y[1 + i]}}, {i, l - 1}];

(soln = Solve[
Join[{eq1, eq2},
eq3 /. {E^(-x h[i_]) :> 1/e[i], E^(x h[i_]) :> e[i]} //
Expand // Simplify],
vars]) // Short // AbsoluteTiming


The solution is quite big and chokes Simplify:

(solnExp = soln /. e[i_] :> E^(x h[i])) // ByteCount // AbsoluteTiming
(*
{15.650537, 729043792}
*)


It's interesting that it takes longer to back substitute, than to solve.

We can simplify piecemeal, working from the bottom up. We can map Simplify to the parts at some level from the bottom (a sort of divide-and-conquer strategy). I use Fold to do it sequentially for levels -6, -8, -10, -12, printing out the timings and bytecounts for each level. The size is quickly reduced.

Depth[solnExp]
(*
38
*)

soln1 = Fold[
(Print[{#[[1]], ByteCount@#[[2, 1]], #[[2, 2]]}]; #[[2, 1]]) &@
AbsoluteTiming @ {Map[Simplify, #1, {#2}], #2} &,
solnExp, -Range[6, 12, 2]]; // AbsoluteTiming

(*   timing    bytecount   level
{18.190976, 465295824,   -6}
{11.706499, 261725136,   -8}
{3.647407,  225441152,  -10}
{17.897693, 166903336,  -12}

{52.811055, Null}
*)

Depth[soln1]
(*
36
*)


We can try to go further up, levels -15, -20, -25, -30, but there are diminishing returns:

soln2 = Fold[
(Print[{#[[1]], ByteCount@#[[2, 1]], #[[2, 2]]}]; #[[2, 1]]) &@
AbsoluteTiming@{Map[Simplify, #1, {#2}], #2} &,
soln1, -Range[15, 30, 5]]; // AbsoluteTiming

(*
{338.471629, 146285856, -15}
hmm…if it ever finishes I'll update.
*)

-