Speed up quad sus tf calculation 10x, by pre-solving the system symbolically

Profiling of precompIFO() before this change:

%timeit gwinc.precompIFO(freq, ifo)

404 ms ± 24 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%prun ifo2 = gwinc.precompIFO(freq, ifo)

      693645 function calls (693211 primitive calls) in 0.833 seconds

Ordered by: internal time

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
 30000    0.362    0.000    0.695    0.000 linalg.py:319(solve)
    10    0.083    0.008    0.778    0.078 suspension.py:39(calc_transfer_functions)
 30008    0.065    0.000    0.109    0.000 linalg.py:141(_commonType)
 30018    0.056    0.000    0.056    0.000 {method 'astype' of 'numpy.ndarray' objects}
 60012    0.046    0.000    0.098    0.000 linalg.py:108(_makearray)
 etc, etc

And after:

%timeit gwinc.precompIFO(freq, ifo)

28.6 ms ± 4.87 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%prun ifo2 = gwinc.precompIFO(freq, ifo)

      3596 function calls (3162 primitive calls) in 0.032 seconds

Ordered by: internal time

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
    74    0.016    0.000    0.016    0.000 coatingthermal.py:743(getCoatRefl)
    10    0.006    0.001    0.006    0.001 suspension.py:30(tst_force_to_tst_displ)
     1    0.005    0.005    0.012    0.012 suspension.py:51(suspQuad)
   137    0.001    0.000    0.001    0.000 {built-in method numpy.core.multiarray.zeros}
     2    0.001    0.000    0.001    0.000 suspension.py:40(top_displ_to_tst_displ)
etc, etc

gwinc/suspension.py

@@ -15,16 +15,34 @@ FIBER_TYPES = [
 ]
 
 
-# quad pendulum equation of motion matrix A =
-# [[k0+k1-m0*w**2,           -k1,             0,          0],
-#  [          -k1, k1+k2-m1*w**2,           -k2,          0],
-#  [            0,           -k2, k2+k3-m2*w**2,        -k3],
-#  [            0,             0,           -k3, k3-m3*w**2]])
-# diagonal elements: mass and restoring forces
-# off-diagonal: coupling to stages above and below
-# want TM equations of motion, so index 4
-# b = [[0], [0], [0], [1]]
-# sympy.linsolve((A, b), x) yields the following two functions
+def generate_symbolic_tfs(stages=4):
+    import sympy as sp
+
+    # construct quad pendulum equation of motion matrix
+    ksyms = sp.numbered_symbols('k')
+    msyms = sp.numbered_symbols('m')
+    w = sp.symbols('w')
+    k = [next(ksyms) for n in range(stages)]
+    m = [next(msyms) for n in range(stages)]
+    A = sp.zeros(stages)
+    for n in range(stages-1):
+        # mass and restoring forces (diagonal elements)
+        A[n, n] = k[n] + k[n+1] - m[n] * w**2
+        # couplings to stages above and below
+        A[n, n+1] = -k[n+1]
+        A[n+1, n] = -k[n+1]
+    # mass and restoring force of bottom stage
+    A[-1, -1] = k[-1] - m[-1] * w**2
+
+    # want TM equations of motion, so index 4
+    b = sp.zeros(stages, 1)
+    b[-1] = 1
+
+    # solve linear system
+    xsyms = sp.numbered_symbols('x')
+    x = [next(xsyms) for n in range(stages)]
+    ans = sp.linsolve((A, b), x)
+    return ans
 
 
 def tst_force_to_tst_displ(k, m, f):
@@ -37,6 +55,7 @@ def tst_force_to_tst_displ(k, m, f):
     X3 = (k2**2*(k0 + k1 - m0*w**2) + (k1**2 - (k0 + k1 - m0*w**2)*(k1 + k2 - m1*w**2))*(k2 + k3 - m2*w**2))/(-k3**2*(k1**2 - (k0 + k1 - m0*w**2)*(k1 + k2 - m1*w**2)) + (k3 - m3*w**2)*(k2**2*(k0 + k1 - m0*w**2) - (-k1**2 + (k0 + k1 - m0*w**2)*(k1 + k2 - m1*w**2))*(k2 + k3 - m2*w**2)))
     return X3
 
+
 def top_displ_to_tst_displ(k, m, f):
     """transfer function for quad pendulum