Added loop stability analysis

idmlin.py

 from scipy import signal
 from scipy import sqrt, linspace, exp, pi, polymul, sinh, cosh, arange
+from scipy import real, array, eye, kron, diag, sparse
+from numpy.linalg import eig
 
 class IDMLin:
     D = 4
@@ -13,27 +15,45 @@ class IDMLin:
     def go(self, v0, re):
         self.v0 = v0
         self.re = re
-        self.f, self.g, self.h, self.ve, self.he = self._fghvh(v0, re)
+        self._fghvh()
 
-    def _fghvh(self, v0, re):
+    def _fghvh(self):
         t = self.t
         a = self.a
         b = self.b
         d = self.D
         s0 = self.s0
+        v0 = self.v0
+        re = self.re
 
         r1 = (1-re**d)
         r32 = r1**(1.5)
         ge = s0 + t * v0 * re
 
-        ve = re * v0
-        he = ge * r1**(-0.5)
+        self.ve = re * v0
+        self.he = ge * r1**(-0.5)
 
-        f = 1*a*r32/ge
-        g = v0 * sqrt(a/b) * re * r1/ge
-        h = a * d * re**(d-1)/v0 + 2*a*t*r1/ge
+        self.f = 1*a*r32/ge
+        self.g = v0 * sqrt(a/b) * re * r1/ge
+        self.h = a * d * re**(d-1)/v0 + 2*a*t*r1/ge
 
-        return f,g,h, ve, he
+    def _mn(self):
+        self.m = array([[0,1],[-self.f, -self.g-self.h]])
+        self.n = array([[0,0],[self.f, self.g]])
+
+    def loop(self, cnt):
+        self._mn()
+        m = self.m
+        n = self.n
+        a = kron(eye(cnt), m) + kron(diag([1]*(cnt-1), -1) + diag([1], cnt-1), n)
+        # a = sparse.kron(sparse.identity(cnt), m) + sparse.kron(sparse.diags(([1]*(cnt-1), [1]), (-1, cnt-1)), n)
+        return a
+
+    def isloopstable(self, cnt):
+        a = self.loop(cnt)
+        w,v = eig(a)
+        # w,v = sparse.linalg.eigs(a)
+        return all(real(w)<1e-10)
 
     def os(self):
         f,g,h = self.f, self.g, self.h
@@ -74,6 +94,12 @@ if __name__ == "__main__":
     idm = IDMLin(t=1.5, a=.3, b=3, s0=2)
     v0 = 30
 
+    re = 0.5
+    idm.go(v0, re)
+    print idm.loop(3)
+    print idm.isloopstable(3)
+    print idm.isloopstable(100)
+    
     res = arange(100)/100.
     idm.go(v0, res)
     plt.plot(res, idm.os())

idmsweep.py

@@ -102,7 +102,6 @@ if __name__ == "__main__":
     sweep(idmgo, rngs, axis=[0,1,-.2,.2])
     '''
 
-    '''
     ### What happens when you vary IDM parameters?
     dots = (':', '--', '-', '.-')
 
@@ -116,8 +115,8 @@ if __name__ == "__main__":
 
     go = linfty(0.5, n=50, tmax=500)
     sweep(go, rngs)
-    '''
 
+    '''
     go = linfty(0.05, n=80, tmax=500, timeplot=True)
     T, Us = go(1.5, 0.3, 3, 2, 30)
 
@@ -125,6 +124,7 @@ if __name__ == "__main__":
     axarr[0].plot(T, Us)
     axarr[1].plot(T, diff(insert(Us, 0, 1, axis=1)))
     plt.show()
+    '''
 
     '''
     go = linfty(0.5, n=500, tmax=2000, dt=1, timeplot=False)