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

import sympy as sp
sp.init_printing()

def generate_lmf(k, lmin_alpha, 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

a,b = generate_lmf(1,0,0,0)

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

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 );

plot_mod1_locus(a,b)

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

k = 1
a,b = generate_lmf(k,k-1,0,k-1)
plot_mod1_locus(a,b)

k = 2
a,b = generate_lmf(k,k-1,0,k-1)
plot_mod1_locus(a,b)

k = 3
a,b = generate_lmf(k,k-1,0,k-1)
plot_mod1_locus(a,b)

k = 4
a,b = generate_lmf(k,k-1,0,k-1)
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]

$\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]

$\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]

generate_lmf(2,0,0,1)   # "full" 2-step explicit

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

$\displaystyle \left( \left[ -5, \ 4, \ 1\right], \ \left[ 2, \ 4, \ 0\right]\right)$

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
    Y = y0 + 0.001*np.random.rand(k)
    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(20): #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.',markersize=3) 
        tj += h
        if abs(Y[-1]) > 100: break
    ax.set_title('h = ' + str(h))

A method that is consistent (p=3) but not stable

# "full" 2-step (uses all available data)
a,b = generate_lmf(2,0,0,1)

def f(t,y): return 0#-y
y0 = 1; t0=0; t1=3
hvals = [0.1,0.05,.025,.0125,.00001]
fig,axes = plt.subplots(1,len(hvals),figsize=(len(hvals)*4,6))
for i,h in enumerate( hvals ):
    use_multistep(axes[i],a,b,f,y0,t0,t1,h)

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

Blows up regardless how small we make $h$, and blows up faster for smaller $h$.

A method that is consistent and stable hence convergent, but not A.S. for any real $h\lambda$

# 2-step midpoint: stable but not A.S. for any real hl
a = [-1,0,1]
b = [0,2,0]

def f(t,y): return -y
y0 = 1; t0=0; t1=3
hvals = [0.1,0.05,.025,.0125]
fig,axes = plt.subplots(1,len(hvals),figsize=(len(hvals)*4,6))
for i,h in enumerate( hvals ):
    use_multistep(axes[i],a,b,f,y0,t0,t1,h)

We see convergence above, $h\rightarrow 0$, but if we extend the timeline, keeping $h$ as in the last plot above ...

ax = plt.subplot(111)
t1 = 10
use_multistep(ax,a,b,f,y0,t0,t1,h)

we see that it blows up eventually.

# full 1-step
generate_lmf(1,0,0,1)

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

$\displaystyle \left( \left[ -1, \ 1\right], \ \left[ \frac{1}{2}, \ \frac{1}{2}\right]\right)$

Trapezoidal method region of absolute stability

def eigenvalue(hlre,hlim):
    hlambda = hlre+hlim*1j # form the complex hlambda from real and imaginary parts
    return abs( (1 + hlambda/2)/(1-hlambda/2) )  # midpoint method

import numpy as np

r = 3  # plot box radius
x = np.arange(-r,r, 0.1)
y = np.arange(-r,r, 0.1)
X, Y = np.meshgrid(x, y)  # grid of points to evaluate eigenvalue

Z = eigenvalue(X,Y)

fig, ax = plt.subplots(figsize=(6,6))
ax.contourf(X,Y,Z,levels=[0,1])
ax.contour(X,Y,Z,levels=[1])
plt.grid()
plt.title('region of absolute stability in $\mathbb{C}$ for method')
plt.show()

The region of absolute stability is the entire left half plane.

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 );

k = 4x`
a,b = generate_lmf(k-1,k,0,k-1) 
plot_mod1_locus(a,b)

$\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]

backward differentiation formulas

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

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

alphas: [-2/11, 9/11, -18/11, 1]
betas : [0, 0, 0, 6/11]

Exterior of red curve is the region of absolute stability: it contains the entire negative real line. Methods are not A-stable, but as good as, for many practical purposes.

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

$\displaystyle \left\{ \alpha_{0} : \frac{3}{25}, \ \alpha_{1} : - \frac{16}{25}, \ \alpha_{2} : \frac{36}{25}, \ \alpha_{3} : - \frac{48}{25}, \ \beta_{4} : \frac{12}{25}\right\}$

alphas: [3/25, -16/25, 36/25, -48/25, 1]
betas : [0, 0, 0, 0, 12/25]