import sympy as sp
sp.init_printing()

A function to generate a linear multistep formula

You can uncomment "display" statements in the code below to see what is going on.

def generate_lmf(lmin_alpha, k, lmin_beta, lmax_beta):

    # create the unknowns
    alpha = [0]*(k+1)
    beta  = [0]*(k+1)
    
    # put sympy symbols in the appropriate slots
    varlist = []
    alpha[k] = 1   # w.o.l.o.g.
    for i in range(lmin_alpha,k): 
        alpha[i] = sp.symbols(f'alpha_{i}')
        varlist.append(alpha[i])
    for i in range(lmin_beta ,lmax_beta+1):
        beta[i]  = sp.symbols(f'beta_{i}')
        varlist.append(beta[i])
    nfree = len(varlist)
    #display(alpha)
    #display(beta)
    #display(nfree)
    #display(varlist)
    
    # build the equations
    c = [sum(alpha)]
    for q in range(1,nfree):
        c.append( sum( [ l**q    *alpha[l]/sp.factorial(q)
                        -l**(q-1)* beta[l]/sp.factorial(q-1) for l in range(k+1) ] ) )
    #for item in c:
    #    display(item)
    #print()
    sol = sp.solve(c,varlist)
    display(sol)
    #print(type(sol))
    #print(type(key) for key in sol)
    
    alpha = [ sol[item] if item in sol else item for item in alpha ]
    beta  = [ sol[item] if item in sol else item for item in beta  ]
    return alpha, beta
    #return [alpha#,beta.subs(sol)
    # can we return a function that implements the corresponding method?
    #def method(t,y,f):

imports and notebook configuration

from matplotlib import rcdefaults
rcdefaults()  # restore default matplotlib rc parameters
%config InlineBackend.figure_format='retina'   
import seaborn as sns  # wrapper for matplotlib that provides prettier styles and more
import matplotlib.pyplot as plt  # use matplotlib functionality directly
%matplotlib inline   
sns.set()

import numpy as np

Complex numbers in Python

$i^2=-1$, for example:

1j**2

(-1+0j)

Plot the locus in the complex $h\lambda$ plane of eigenvalue of modulus 1

I also draw the locus of eigenvalues of modulus 1-epsilon (in faint green) to suggest which side of the boundary is the |z| < 1 side.

import numpy as np
def plot_mod1_locus(a,b):
    print('alphas:',a)
    print('betas :',b)
    a = np.array(a,dtype=float)
    b = np.array(b,dtype=float)
    k = len(a)-1

    plt.subplot(111,aspect=1)
    z = np.exp( (0+1j)*np.linspace(0,2*np.pi,200) )
    #print(z)
    plt.axhline(0,color='gray')
    plt.axvline(0,color='gray')
    def get_hlambda(z): return np.sum([a[l]*z**l for l in range(k+1)],axis=0)/np.sum([b[l]*z**l for l in range(k+1)],axis=0)
    hlambda = get_hlambda(z)
    plt.plot( np.real(hlambda),np.imag(hlambda),'r' );

    hlambda = get_hlambda(0.965*z)  # draw the locus of eigenvalues of modulus 1-epsilon to indicate where |z| < 1
    plt.plot( np.real(hlambda),np.imag(hlambda),'g',alpha=0.45 );

Example: Adams-Bashforth 4-step

a,b = generate_lmf( 3,4,0,3) # Adams-Bashforth 4 
plot_mod1_locus(a,b)
plt.title('Adams-Bashforth 4'+' stability region');

$\displaystyle \left\{ \alpha_{3} : -1, \ \beta_{0} : - \frac{3}{8}, \ \beta_{1} : \frac{37}{24}, \ \beta_{2} : - \frac{59}{24}, \ \beta_{3} : \frac{55}{24}\right\}$

alphas: [0, 0, 0, -1, 1]
betas : [-3/8, 37/24, -59/24, 55/24, 0]

The semi-disc to the left of the imaginary axis is the stability region where all 4 eigenvalues have modulus < 1.

Outside of that, there is at least one eigenvalue with modulus > 1, and in the "ears" there are two of them.

A function to apply the multistep method characterized by the given alphas and betas

to the IVP y' = f(t,y), y(t0)=y0, up to t=t1 using step-size h

def use_multistep(ax,alpha,beta,f,y0,t0,t1,h):  
    #print('alphas:' + str(alpha))
    #print('betas: ' + str(beta))
    assert( beta[-1] == 0 ) # this works for explicit methods only: implicit requires more work
    k = len(alpha)-1
    Y = np.array([y0 + l*h*f(t0,y0) for l in range(k)])  # crude Euler initialization of state vector
    F = np.array([f(t0+l*h,Y[l])    for l in range(k)])
    #print(Y)
    #print(F)
    ax.plot( [t0+l*h for l in range(k)], Y, '.')
    tj = t0
    ax.axhline(0,color='gray',alpha=0.5)
    for j in range(int((t1-t0)/h)):
        y_j_plus_k = - np.sum( [alpha[l]*Y[l] - beta[l]*h*F[l] for l in range(k)] )
        Y[:k-1] = Y[1:]
        F[:k-1] = F[1:]
        Y[k-1] = y_j_plus_k
        F[k-1] = f(tj+k*h,Y[k-1]) # our one function evaluation per step
        ax.plot(tj+k*h,Y[k-1],'r.') 
        tj += h
        if abs(Y[-1]) > 100: break
    ax.set_title('h = ' + str(h))

Euler appied to y'=-y inside and outside the stability region

a,b = generate_lmf( 0,1,0,0) # Euler
plot_mod1_locus(a,b)
plt.title('Euler'+' stability region');

$\displaystyle \left\{ \alpha_{0} : -1, \ \beta_{0} : 1\right\}$

alphas: [-1, 1]
betas : [1, 0]

hvals = [0.1,0.5,1.5,2.5]
fig,axes = plt.subplots(1,len(hvals),figsize=(len(hvals)*3,3))
for i,h in enumerate( hvals ):
    use_multistep(axes[i],a,b,f,y0,t0,t1,h)

We see the solution gets increasingly poor as we approach the boundary, and blows up once we are no longer inside it.

Adams-Bashforth 4-step appied to y'=-y inside and outside the stability region

a,b = generate_lmf( 3,4,0,3) # Adams-Bashforth 4 
plot_mod1_locus(a,b)
plt.title('Adams-Bashforth 4'+' stability region');

$\displaystyle \left\{ \alpha_{3} : -1, \ \beta_{0} : - \frac{3}{8}, \ \beta_{1} : \frac{37}{24}, \ \beta_{2} : - \frac{59}{24}, \ \beta_{3} : \frac{55}{24}\right\}$

$\displaystyle \left[ 0, \ 0, \ 0, \ -1, \ 1\right]$

$\displaystyle \left[ - \frac{3}{8}, \ \frac{37}{24}, \ - \frac{59}{24}, \ \frac{55}{24}, \ 0\right]$

Note that the interval of stability on the real axis is much narrower than for Euler's method (though the method has higher order, 4).

def f(t,y): 
    return -y

hvals = [.1,.25,.30,.31]
fig,axes = plt.subplots(1,len(hvals),figsize=(len(hvals)*3,3))
for i,h in enumerate( hvals ):
    use_multistep(axes[i],a,b,f,y0,t0,t1,h)

Again, the solution blows up when we venture outside of the stability region.

Stability regions for Adams-Bashforth family k = 1,2,3,4,...

for k in [1,2,3,4]:
    a,b = generate_lmf( k-1,k,0,k-1) # Adams-Bashforth k-step
    plot_mod1_locus(a,b)

$\displaystyle \left\{ \alpha_{0} : -1, \ \beta_{0} : 1\right\}$

alphas: [-1, 1]
betas : [1, 0]

$\displaystyle \left\{ \alpha_{1} : -1, \ \beta_{0} : - \frac{1}{2}, \ \beta_{1} : \frac{3}{2}\right\}$

alphas: [0, -1, 1]
betas : [-1/2, 3/2, 0]

<ipython-input-141-d0dc759f2cd1>:9: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  plt.subplot(111,aspect=1)

$\displaystyle \left\{ \alpha_{2} : -1, \ \beta_{0} : \frac{5}{12}, \ \beta_{1} : - \frac{4}{3}, \ \beta_{2} : \frac{23}{12}\right\}$

alphas: [0, 0, -1, 1]
betas : [5/12, -4/3, 23/12, 0]

<ipython-input-141-d0dc759f2cd1>:9: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  plt.subplot(111,aspect=1)

$\displaystyle \left\{ \alpha_{3} : -1, \ \beta_{0} : - \frac{3}{8}, \ \beta_{1} : \frac{37}{24}, \ \beta_{2} : - \frac{59}{24}, \ \beta_{3} : \frac{55}{24}\right\}$

alphas: [0, 0, 0, -1, 1]
betas : [-3/8, 37/24, -59/24, 55/24, 0]

<ipython-input-141-d0dc759f2cd1>:9: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  plt.subplot(111,aspect=1)