Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

enter image description here

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

share|improve this question
    
I have deleted Sjoerd's answer at his request because he stated that it was incorrect. –  Mr.Wizard Dec 31 '13 at 17:36
add comment

1 Answer

up vote 4 down vote accepted

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

Mathematica graphics

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.
*)
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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