Stability of Euler's method

# first get resources306 module
import requests
r = requests.get('https://raw.githubusercontent.com/UBmath/306/master/resources306.py')
with open('resources306.py','w') as f: f.write(r.text)

from resources306 import *
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()

def euler(f,t0,y0,h,m):  # written as an RK method
    
    b21 = 1
    
    t = t0
    y = y0

    ya = np.empty(m+1)
    ya[0] = y
    ta = np.linspace(t0,t0+m*h,m+1)
     
    for k in range(m):
        t = t0 + k*h
        K1 = f(t,y)  
        y1 = y + h*b21*K1
        y = y1
        ya[k+1] = y
        
    return ta,ya

def f(t,y): return -y

plt.figure(figsize=(20,10))
slopefieldplot(f,0,15,-2,2,.2)
for h in [2.1,1.8,.9,.2]:
    ta,ya = euler(f,0,1.0,h,int(15/h))
    plt.plot(ta,ya,'o-',clip_on=False,label='h='+str(h))
plt.legend()
plt.xlabel('t'); plt.ylabel('y');
plt.title('Absolute stability: Euler applied to $y^\prime = -y$')
plt.savefig('day08_1_euler_stability.png',bbox_inches='tight')

Stability region by countour plot of $|\mu|$

from resources306 import *
%config InlineBackend.figure_format='retina'   

def mu(hlambda):
    return 'euler', 1 + hlambda

r = 3  # plot box radius
re = np.arange(-r,r, 0.1)
im = np.arange(-r,r, 0.1)
RE, IM = np.meshgrid(re, im)  # grid of points to evaluate eigenvalue
hlambda = RE + 1j*IM

method,MU = mu(hlambda)

fig, ax = plt.subplots(figsize=(6,6))
ax.contourf(RE,IM,np.abs(MU),levels=[0,1])
ax.contour (RE,IM,np.abs(MU),levels=[1])
plt.axvline(0)
plt.axhline(0)
plt.grid(True)
plt.title('region of absolute stability in $\mathbb{C}$ for ' + method + ' method')
plt.text(0.8,1.5,'$h\lambda \in \mathbb{C}$',fontsize=20)
plt.savefig(f'day08_2_stability_region_{method}.png',bbox_inches='tight')

from resources306 import *
%config InlineBackend.figure_format='retina'   

def mu(hlambda):
    return 'midpoint', 1 + hlambda + hlambda**2/2

r = 3  # plot box radius
re = np.arange(-r,r, 0.1)
im = np.arange(-r,r, 0.1)
RE, IM = np.meshgrid(re, im)  # grid of points to evaluate eigenvalue
hlambda = RE + 1j*IM

method,MU = mu(hlambda)

fig, ax = plt.subplots(figsize=(6,6))
ax.contourf(RE,IM,np.abs(MU),levels=[0,1])
ax.contour (RE,IM,np.abs(MU),levels=[1])
plt.axvline(0)
plt.axhline(0)
plt.grid(True)
plt.title('region of absolute stability in $\mathbb{C}$ for ' + method + ' method')
plt.text(0.8,1.5,'$h\lambda \in \mathbb{C}$',fontsize=20)
plt.savefig(f'day08_2_stability_region_{method}.png',bbox_inches='tight')

Stability region for RK4

import sympy as sp
from IPython.display import HTML

h,lamda = sp.symbols('h,lambda')
def f(t,y):
    return lamda*y

def rk4(f,t,y,h):

    b21 = sp.Rational(1,2)
    b31 = 0
    b32 = sp.Rational(1,2)
    b41 = 0
    b42 = 0
    b43 = 1
    b51 = sp.Rational(1,6)
    b52 = sp.Rational(1,3)
    b53 = sp.Rational(1,3)
    b54 = sp.Rational(1,6)
        
    K1 = f(t,y)  
    y1 = y + h*b21*K1
    K2 = f(t+h*b21,y1)
    y2 = y + h*(b31*K1 + b32*K2)
    K3 = f(t+h*(b31    + b32   ), y2)
    y3 = y + h*(b41*K1 + b42*K2 + b43*K3)
    K4 = f(t+h*(b31    + b32    + b43   ), y3)
    y4 = y + h*(b51*K1 + b52*K2 + b53*K3 + b54*K4)
    y = y4
            
    return y

mu = rk4(f,0,1,h).expand()
display(HTML('$\mu =$'),mu)

r = 3.1  # plot box radius
mu = sp.lambdify(lamda,rk4(f,0,1,1).expand())

re = np.arange(-r,r, 0.1)
im = np.arange(-r,r, 0.1)
RE, IM = np.meshgrid(re, im)  # grid of points to evaluate eigenvalue
hlambda = RE + 1j*IM

method = 'RK4'
MU = mu(hlambda)

fig, ax = plt.subplots(figsize=(6,6))
ax.contourf(RE,IM,np.abs(MU),levels=[0,1])
ax.contour (RE,IM,np.abs(MU),levels=[1])
plt.axvline(0)
plt.axhline(0)
plt.grid(True)
plt.title('region of absolute stability in $\mathbb{C}$ for ' + method + ' method')
plt.text(0.8,1.5,'$h\lambda \in \mathbb{C}$',fontsize=20)
plt.savefig(f'day08_2_stability_region_{method}.png',bbox_inches='tight')

$\mu =$

$\displaystyle \frac{h^{4} \lambda^{4}}{24} + \frac{h^{3} \lambda^{3}}{6} + \frac{h^{2} \lambda^{2}}{2} + h \lambda + 1$

RK4 applied to $y'=-y$ for various $h$

def use_rk4(f,t0,y0,h,m):

    b21 = 1/2
    b31 = 0
    b32 = 1/2
    b41 = 0
    b42 = 0
    b43 = 1
    b51 = 1/6
    b52 = 1/3
    b53 = 1/3
    b54 = 1/6
        
    t = t0
    y = y0

    ya = np.empty(m+1)
    ya[0] = y
    ta = np.linspace(t0,t0+m*h,m+1)
     
    for k in range(m):
        t = t0 + k*h
        K1 = f(t,y)  
        y1 = y + h*b21*K1
        K2 = f(t+h*b21,y1)
        y2 = y + h*(b31*K1 + b32*K2)
        K3 = f(t+h*(b31    + b32   ), y2)
        y3 = y + h*(b41*K1 + b42*K2 + b43*K3)
        K4 = f(t+h*(b31    + b32    + b43   ), y3)
        y4 = y + h*(b51*K1 + b52*K2 + b53*K3 + b54*K4)
        y = y4
        ya[k+1] = y
        
    return ta,ya

def f(t,y):
    return -y

for h in [.1,1,2.7,2.785293563405282,2.8]:
    ta,ya = use_rk4(f,0,1.0,h,int(10/h))
    plt.plot(ta,ya,'-o')
ta = np.linspace(0,10,100)
plt.plot(ta,np.exp(-ta),'k'); plt.xlabel('t'); plt.ylabel('y',rotation=0)
plt.grid()
plt.savefig(f'day08_3_RK4_applied_various_h.png',bbox_inches='tight');

Comparison with exact solution

from resources306 import *
h,lamda = sp.symbols('h,lambda')
def f(t,y):
    return lamda*y
hl = sp.symbols('hlambda')
expressionplot( 1,hl,-3,1)
expressionplot( 0,hl,-3,1)
expressionplot(-1,hl,-3,1)

expressionplot(1 + hl ,hl,-3,1,label='Euler')
expressionplot(1 + hl + hl**2/2  ,hl,-3,1,label='midpoint')
expressionplot(1 + hl + hl**2/2 + hl**3/6 + hl**4/24 ,hl,-3,1,label='RK4')
expressionplot(sp.exp(hl),hl,-3,1,color='k',alpha=.8,label='exact')
plt.axvline(0,color='k')
plt.legend()
plt.savefig(f'day08_4_comparison_with_exact.png',bbox_inches='tight');

Generate linear multistep formula using specified data

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

import pandas as pd
def nice(lmf): display(pd.DataFrame(lmf,index=('alpha','beta')))

# Euler
euler = generate_lmf(1,0,0,0)
nice(euler)

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

	0	1
alpha	-1	1
beta	1	0

# Adams-Bashforth 4th order
ab4 = generate_lmf(4,3,0,3)
nice(ab4)

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

	0	1	2	3	4
alpha	0	0	0	-1	1
beta	-3/8	37/24	-59/24	55/24	0

compare the above with Wikipedia on Linear Multistep Methods

[5]*3

$\displaystyle \left[ 5, \ 5, \ 5\right]$