In [137]:
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.

In [138]:
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¶

In [139]:
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:

In [140]:
1j**2
Out[140]:
(-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.

In [141]:
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¶

In [142]:
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

In [143]:
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¶

In [144]:
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]
In [136]:
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¶

In [111]:
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).

In [110]:
def f(t,y):
return -y
In [109]:
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,...¶

In [146]:
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)
In [ ]: