Formatting, bugfixes, and a feature.

Several PEP8 formatting corrections.
Remove calls to plt.hold.
Add 3rd-order additive LMMs from ARW paper.
Miscellaneous bug fixes.

.travis.yml

@@ -15,3 +15,10 @@ script: coverage run --source=nodepy tests.py
 after_success:
   - coveralls
   - codecov
+
+deploy:
+  provider: pypi
+  user: ketch
+  password: ...
+  on:
+    tags: true

nodepy/ivp.py

@@ -85,21 +85,24 @@ def load_ivp(ivpname='All'):
 
     name = 'test'
     description = r"Dahlquist's test problem; $f(y) = \lambda y$"
+    u0 = 1.
     rhs = lambda t,u: u
-    exact = lambda t : ivp.u0*np.exp(t)
-    ivps[name] = IVP(f=rhs, u0=1., T=5., desc=description, exact=exact, name=name)
+    exact = lambda t: u0*np.exp(t)
+    ivps[name] = IVP(f=rhs, u0=u0, T=5., desc=description, exact=exact, name=name)
 
     name = 'zoltan'
+    u0 = 1.
     rhs = lambda t,u: -100*abs(u)
-    exact = lambda t : ivp.u0*np.exp(-100*t)
+    exact = lambda t: u0*np.exp(-100*t)
     description = 'The linear scalar test problem with abs value'
-    ivps[name] = IVP(f=rhs, u0=1., T=10., desc=description, exact=exact, name=name)
+    ivps[name] = IVP(f=rhs, u0=u0, T=10., desc=description, exact=exact, name=name)
 
     name = 'nlsin'
+    u0 = 1.
     rhs = lambda t,u: 4.*u*float(np.sin(t))**3*np.cos(t)
-    exact = lambda t: ivp.u0*np.exp((np.sin(t))**4)
+    exact = lambda t: u0*np.exp((np.sin(t))**4)
     description = 'A simple nonlinear scalar problem'
-    ivps[name] = IVP(f=rhs, u0=1., T=5., desc=description, exact=exact, name=name)
+    ivps[name] = IVP(f=rhs, u0=u0, T=5., desc=description, exact=exact, name=name)
 
     name = 'ode1'
     rhs = lambda t,u: 4.*t*np.sqrt(u)

nodepy/linear_multistep_method.py

@@ -36,9 +36,7 @@
 from six.moves import range
 
 
-#=====================================================
 class LinearMultistepMethod(GeneralLinearMethod):
-#=====================================================
     r"""
         Implementation of linear multistep methods in the form:
 
@@ -162,11 +160,11 @@ def latex(self):
         for i in range(k+1):
             subscript = latex(n+k-i)
             if self.alpha[k-i] != 0:
-                alpha_terms.append(printable(self.alpha[k-i],return_one=False)
-                                   + ' y_{'+subscript+'}')
+                alpha_terms.append(printable(self.alpha[k-i],return_one=False) +
+                                   ' y_{'+subscript+'}')
             if self.beta[k-i] != 0:
-                beta_terms.append(printable(self.beta[k-i],return_one=False)
-                                  +  'h f(y_{'+subscript+'})')
+                beta_terms.append(printable(self.beta[k-i],return_one=False) +
+                                  'h f(y_{'+subscript+'})')
         lhs = ' + '.join(alpha_terms)
         rhs = ' + '.join(beta_terms)
         s = r'\begin{align}'+ ' = '.join([lhs,rhs]) + r'\end{align}'
@@ -201,7 +199,6 @@ def ssp_coefficient(self):
         return min([-self.alpha[j]/self.beta[j]
                     for j in range(len(self.alpha)-1) if self.beta[j]!=0])
 
-
     def plot_stability_region(self,N=100,bounds=None,color='r',filled=True, alpha=1.,
                               to_file=False,longtitle=False):
         r"""
@@ -233,7 +230,7 @@ def plot_stability_region(self,N=100,bounds=None,color='r',filled=True, alpha=1.
         """
         import matplotlib.pyplot as plt
         rho, sigma = self.__num__().characteristic_polynomials()
-        mag = lambda z : _root_condition(rho-z*sigma)
+        mag = lambda z: _root_condition(rho-z*sigma)
         vmag = np.vectorize(mag)
 
         z = self._boundary_locus()
@@ -253,7 +250,6 @@ def plot_stability_region(self,N=100,bounds=None,color='r',filled=True, alpha=1.
 
         R=1.5-vmag(Z)
 
-
         if filled:
             plt.contourf(X,Y,R,[0,1],colors=color,alpha=alpha)
         else:
@@ -266,12 +262,10 @@ def plot_stability_region(self,N=100,bounds=None,color='r',filled=True, alpha=1.
         else:
             plt.setp(ax,title='Stability region')
 
-        plt.hold(True)
         plt.plot([0,0],[bounds[2],bounds[3]],'--k',linewidth=2)
         plt.plot([bounds[0],bounds[1]],[0,0],'--k',linewidth=2)
         plt.plot(np.real(z),np.imag(z),color='k',linewidth=3)
         plt.axis(bounds)
-        plt.hold(False)
 
         if to_file:
             plt.savefig(to_file, transparent=True, bbox_inches='tight', pad_inches=0.3)
@@ -298,15 +292,12 @@ def plot_boundary_locus(self,N=1000):
         plt.figure()
         plt.plot(np.real(z),np.imag(z),color='k',linewidth=3)
         plt.axis('image')
-        plt.hold(True)
         bounds = plt.axis()
         plt.plot([0,0],[bounds[2],bounds[3]],'--k',linewidth=2)
         plt.plot([bounds[0],bounds[1]],[0,0],'--k',linewidth=2)
         plt.title('Boundary locus for '+self.name)
-        plt.hold(False)
         plt.draw()
 
-
     def _boundary_locus(self, N=1000):
         r"""Compute the boundary locus, which is
             given by the set of points
@@ -368,9 +359,7 @@ def __len__(self):
         return len(self.alpha)-1
 
 
-#=====================================================
 class AdditiveLinearMultistepMethod(GeneralLinearMethod):
-#=====================================================
     r"""
         Method for solving equations of the form
 
@@ -405,7 +394,7 @@ def __num__(self):
         numself = copy.deepcopy(self)
         if self.alpha.dtype == object:
             numself.alpha = np.array(self.alpha, dtype=np.float64)
-            numself.beta  = np.array(self.beta,  dtype=np.float64)
+            numself.beta  = np.array(self.beta, dtype=np.float64)
             numself.gamma = np.array(self.gamma, dtype=np.float64)
         return numself
 
@@ -426,8 +415,8 @@ def order(self,tol=1.e-10):
 
         return min(orders)
 
-
-    def plot_imex_stability_region(self,N=100,color='r',filled=True, alpha=1.,fignum=None):
+    def plot_imex_stability_region(self,both_real=False,N=100,color='r',filled=True,
+                                   alpha=1.,fignum=None,bounds=[-10, 1, -5, 5]):
         r"""
             **Input**: (all optional)
                 - N       -- Number of gridpoints to use in each direction
@@ -438,20 +427,20 @@ def plot_imex_stability_region(self,N=100,color='r',filled=True, alpha=1.,fignum
         import matplotlib.pyplot as plt
         rho, sigma1 = self.method1.__num__().characteristic_polynomials()
         rho, sigma2 = self.method2.__num__().characteristic_polynomials()
-        mag = lambda a, b : _max_root(rho - a*sigma1 - 1j*b*sigma2)
+        if both_real:
+            mag = lambda a, b: _max_root(rho - a*sigma1 - b*sigma2)
+        else:
+            mag = lambda a, b: _max_root(rho - a*sigma1 - 1j*b*sigma2)
         vmag = np.vectorize(mag)
-        bounds = [-10, 1, -5, 5]
 
         y = np.linspace(bounds[2],bounds[3],N)
         Y = np.tile(y[:,np.newaxis],(1,N))
         x = np.linspace(bounds[0],bounds[1],N)
         X = np.tile(x,(N,1))
-        Z = X+Y*1j
 
         R = vmag(X,Y)
 
         h = plt.figure(fignum)
-        plt.hold(True)
         if filled:
             plt.contourf(X,Y,R,[0,1],colors=color,alpha=alpha)
         else:
@@ -461,7 +450,6 @@ def plot_imex_stability_region(self,N=100,color='r',filled=True, alpha=1.,fignum
         plt.plot([0,0],[bounds[2],bounds[3]],'--k',linewidth=2)
         plt.plot([bounds[0],bounds[1]],[0,0],'--k',linewidth=2)
         plt.axis(bounds)
-        plt.hold(False)
         return h
 
     def stiff_damping_factor(self,epsilon=1.e-7):
@@ -472,7 +460,7 @@ def stiff_damping_factor(self,epsilon=1.e-7):
         """
         rho, sigma1 = self.method1.__num__().characteristic_polynomials()
         rho, sigma2 = self.method2.__num__().characteristic_polynomials()
-        mag = lambda a, b : _max_root(rho - a*sigma1 - 1j*b*sigma2)
+        mag = lambda a, b: _max_root(rho - a*sigma1 - 1j*b*sigma2)
 
         f=[]
         z=-1.
@@ -487,7 +475,7 @@ def stiff_damping_factor(self,epsilon=1.e-7):
                 raise Exception('Unable to compute stiff damping factor: slow convergence')
 
 
-#======================================================
+# ======================================================
 def _max_root(p):
     return max(np.abs(p.r))
 
@@ -514,9 +502,9 @@ def _root_condition(p,tol=1.e-13):
     return True
 
 
-#======================================================
+# ======================================================
 # Families of multistep methods
-#======================================================
+# ======================================================
 
 def Adams_Bashforth(k):
     r"""
@@ -775,7 +763,6 @@ def sand_cc(s):
     zero = sympy.Rational(0)
 
     k = 2**s + 1
-    p = 2*(s+1)
 
     Jn = [k,k-1]
     for i in range(1,s+1):
@@ -797,7 +784,7 @@ def sand_cc(s):
 
 def arw2(gam,c):
     r"""Returns the second order IMEX additive multistep method based on the
-    parametrization in Section 3.2 of Ascher, Ruuth, & Spiteri.  The parameters
+    parametrization in Section 3.2 of Ascher, Ruuth, & Whetton.  The parameters
     are gam and c.  Known methods are obtained with the following values:
 
         (1/2,0):   CNAB
@@ -821,10 +808,29 @@ def arw2(gam,c):
     """
     half = sympy.Rational(1,2)
     alpha = snp.array([gam-half,-2*gam,gam+half])
-    beta  = snp.array([c/2,1-gam-c,gam+c/2]) #implicit part
-    gamma = snp.array([-gam,gam+1,0]) #explicit part
+    beta  = snp.array([c/2,1-gam-c,gam+c/2])  # implicit part
+    gamma = snp.array([-gam,gam+1,0])  # explicit part
     return AdditiveLinearMultistepMethod(alpha,beta,gamma,'ARW2('+str(gam)+','+str(c)+')')
 
+def arw3(gam,theta,c):
+    r"""Returns the third order IMEX additive multistep method based on the
+    parametrization in Section 3.3 of Ascher, Ruuth, & Whetton.  The parameters
+    are gamma, theta, and c.  Known methods are obtained with the following values:
+
+        (1,0,0):     SBDF3
+
+    Note that there is one sign error in the ARW paper; it is corrected here.
+    """
+    half = sympy.Rational(1,2)
+    third = sympy.Rational(1,3)
+    alpha = snp.array([-half*gam**2+third/2, 3*half*gam**2+gam-1, -3*half*gam**2-2*gam+half-theta,
+                       half*gam**2+gam+third+theta])
+    beta  = snp.array([5*half/6*theta-c, (gam**2-gam)*half+3*c-4*theta*third, 1-gam**2-3*c+23*theta*third/4,
+                       (gam**2+gam)*half+c])
+    gamma = snp.array([(gam**2+gam)*half+5*theta*third/4, -gam**2-2*gam-4*theta*third,
+                       (gam**2+3*gam)*half+1+23*theta*third/4,0])
+    return AdditiveLinearMultistepMethod(alpha,beta,gamma,'ARW3('+str(gam)+','+str(theta)+','+str(c)+')')
+
 
 def loadLMM(which='All'):
     """
@@ -838,17 +844,17 @@ def loadLMM(which='All'):
         True
     """
     LM={}
-    #================================================
+    # ================================================
     alpha = snp.array([-12,75,-200,300,-300,137])/sympy.Rational(137,1)
     beta = snp.array([60,-300,600,-600,300,0])/sympy.Rational(137,1)
     LM['eBDF5'] = LinearMultistepMethod(alpha,beta,'eBDF 5')
-    #================================================
+    # ================================================
     theta = sympy.Rational(1,2)
     alpha = snp.array([-1,1])
     beta = snp.array([1-theta,theta])
     gamma = snp.array([1,0])
     LM['ET112'] = AdditiveLinearMultistepMethod(alpha,beta,gamma,'Euler-Theta')
-    #================================================
+    # ================================================
     if which=='All':
         return LM
     else:

nodepy/ode_solver.py

@@ -101,6 +101,7 @@ def __call__(self,ivp,t0=0,N=5000,dt=None,errtol=None,controllertype='P',
             if errtol is None:      # Fixed-timestep mode
                 if dt is None: dt = (t_final-t0)/float(N)
                 dt_standard = dt + 0
+                max_steps = max(max_steps, int(round(2*(t_final-t0)/dt)))
                 for istep in range(max_steps):
 
                     if t_current+dt >= next_output_time:

nodepy/runge_kutta_method.py

@@ -1670,7 +1670,6 @@ def internal_stability_plot(self,bounds=None,N=200,use_butcher=False,formula='lt
         fig = plt.figure()
         CS = plt.contour(X,Y,th_max,colors='k',levels=levels)
         plt.clabel(CS, fmt='%d', colors='k')#,manual=True)
-        plt.hold(True)
 
         p,q=self.__num__().stability_function(mode='float')
         stability_function.plot_stability_region(p,q,N,color='k',filled=False,bounds=bounds,
@@ -2092,7 +2091,6 @@ def plot_stability_region(self,N=200,color='r',filled=True,bounds=None,
         fig = stability_function.plot_stability_region(p,q,N,color,filled,
                         bounds,plotroots,alpha,scalefac,fignum)
 
-        plt.hold(True)
         p,q = self.embedded_method.__num__().stability_function(mode='float')
         stability_function.plot_stability_region(p,q,N,color='k',filled=False,bounds=bounds,
                 plotroots=plotroots,alpha=alpha,scalefac=scalefac,fignum=fig.number)
@@ -3634,7 +3632,6 @@ def linearly_stable_step_size(rk, L, acc=1.e-7, tol=1.e-14, plot=1):
     h=bisect(0,hmax,acc,tol,_is_linearly_stable, params=(p,q,lamda))
     if plot:
         rk.plot_stability_region()
-        plt.hold(True)
         plt.plot(np.real(h*lamda), np.imag(h*lamda),'o')
     return h
 

nodepy/stability_function.py

@@ -27,10 +27,11 @@ def E_polynomial(p, q=None):
 
     for i in range(1,len(p.c)+1):
         piy.coeffs[-i] = (sympy.I)**(i-1)*piy.coeffs[-i]
-        qiy.coeffs[-i] = (sympy.I)**(i-1)*qiy.coeffs[-i]
-
         pmiy.coeffs[-i] = (-sympy.I)**(i-1)*pmiy.coeffs[-i]
+
+    for i in range(1,len(q.c)+1):
         qmiy.coeffs[-i] = (-sympy.I)**(i-1)*qmiy.coeffs[-i]
+        qiy.coeffs[-i] = (sympy.I)**(i-1)*qiy.coeffs[-i]
 
     Epoly = qiy*qmiy - piy*pmiy
     return Epoly
@@ -208,7 +209,6 @@ def plot_stability_region(p,q,N=200,color='r',filled=True,bounds=None,
 
     # Plot
     h = plt.figure(fignum)
-    plt.hold(True)
     if filled:
         plt.contourf(X,Y,R,[0,1],colors=color,alpha=alpha)
     else:
@@ -218,7 +218,6 @@ def plot_stability_region(p,q,N=200,color='r',filled=True,bounds=None,
     plt.plot([0,0],[bounds[2],bounds[3]],'--k',linewidth=2)
     plt.plot([bounds[0],bounds[1]],[0,0],'--k',linewidth=2)
     plt.axis('Image')
-    plt.hold(False)
     return h
 
 
@@ -254,12 +253,10 @@ def plot_order_star(p,q,N=200,bounds=[-5,5,-5,5], plotroots=False,
 
     h = plt.figure(fignum)
     plt.contourf(X,Y,R,[0,1,1.e299],colors=color)
-    plt.hold(True)
     if plotroots: plt.plot(np.real(p.r),np.imag(p.r),'ok')
     plt.plot([0,0],[bounds[2],bounds[3]],'--k')
     plt.plot([bounds[0],bounds[1]],[0,0],'--k')
     plt.axis('Image')
-    plt.hold(False)
     return h