Runge Kutta methods and stability

import numpy as np
from scipy.optimize import root
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import plotly.io as pio
pio.templates.default = "seaborn"

Curtiss and Hirschfelder problem¶

We consider the following problem :

\left\{ \begin{aligned} {\mathrm d}_t y(t) & = k \, \big(cos(t) - y(t)) \big) \quad \text{with} \quad k > 1\\ y(0) & = y_0 \end{aligned} \right.

and in the following, we will assume $t_0=0$ .

class curtiss_model:

    def __init__(self, k):
        self.k = k

    def fcn(self, t, u) :
        k = self.k
        u_dot = k * (np.cos(t) - u)
        return u_dot

    def sol(self, uini, t0, t):
        k = self.k

        c0 = (uini - (k/(k*k + 1)) * (k*np.cos(t0) + np.sin(t0))) * np.exp(k*t0)
        u = (k/(k*k + 1)) * (k*np.cos(t) + np.sin(t)) +  c0 * np.exp(-k*t)
        return u

Runge Kutta methods¶

def plot_sol_and_error(uini, tini, tend, k, order):
    
    if order==1:
        method_str = "Forward Euler"
        method = forward_euler
        stability_disk = stabiliy_disk_rk1
    elif order==2:
        method_str = "RK2"
        method = rk2
        stability_disk = stabiliy_disk_rk2
    elif order==3:
        method_str = "RK3"
        method = rk3
        stability_disk = stabiliy_disk_rk3
    elif order==4:
        method_str = "RK4"
        method = rk4
        stability_disk = stabiliy_disk_rk4
    else:
        print("Order not implemented")
        exit()
        
    cm = curtiss_model(k)
    fcn = cm.fcn

    texa = np.linspace(tini, tend, 500)
    uexa = cm.sol(uini, tini, texa)
    x, y, disk = stability_disk()
    
    specs=[[{"colspan": 2}, None], [{}, {}]]
    titles=("Solution","Global error", "Stability domain")
    fig = make_subplots(rows=2, cols=2, specs=specs, vertical_spacing=0.15, subplot_titles=titles)
    
    fig.add_trace(go.Scattergl(x=texa, y=uexa, name='Exact solution', legendgroup='1'), row=1, col=1)
    fig.add_trace(go.Contour(x=x, y=y, z=disk, showscale=False, colorscale='blues'), row=2, col=2)

    nt = np.geomspace(25, 200, num=15, dtype=int)
    for nt_i in nt:
        dt = (tend-tini) / (nt_i-1)
        sol = method(tini, tend, nt_i, uini, fcn)
        uexa = cm.sol(uini, tini, sol.t)
        fig.add_trace(go.Scattergl(visible=False, x=sol.t, y=sol.y[0], mode='markers+lines', line_dash='dot', name=method_str,
                                   marker=dict(size=5, symbol='x'), legendgroup='1'), row=1, col=1)
        fig.add_trace(go.Scattergl(visible=False, x=sol.t[1:], y=np.abs(sol.y[0]-uexa)[1:], mode='markers',  showlegend=False,
                                   marker=dict(size=5, symbol='x'),line_color='rgb(221,132,82)', legendgroup='2'), row=2, col=1)
        fig.add_trace(go.Scatter(visible=False, x=[-k*dt], y=[0], legendgroup='3', name='Eigen value',
                                 line_color='rgb(196,78,82)', mode='markers', marker=dict(size=8, symbol='x')), row=2, col=2)

    i_beg = nt.size//2
    fig.data[3*i_beg+2].visible = True
    fig.data[3*i_beg+3].visible = True
    fig.data[3*i_beg+4].visible = True

    # Create and add slider
    steps = []
    for i, nt_i in enumerate(nt):
        args = args=[{"visible": [el==0 or el==1 or el==3*i+2 or el==3*i+3 or el==3*i+4 for el in range(len(fig.data))]}]
        step = dict(method="update", label = f"{nt_i}", args=args)
        steps.append(step)
    sliders = [dict(active=i_beg, currentvalue={"prefix": "nt : "}, steps=steps)]

    legend=dict(x=0.78,y=0.98, bgcolor='rgba(0,0,0,0)', tracegroupgap=270)
    fig.update_xaxes(title='t', col=1)
    fig.update_yaxes(row=1, col=1, title='u')
    fig.update_yaxes(row=2, col=1, title='|error|')
    # fig.update_yaxes(range=[-2,2], row=2, col=2)
    fig.update_layout(sliders=sliders, title=method_str, height=800, legend=legend, modebar=dict(orientation='v'))
    fig['layout']['sliders'][0]['pad']=dict(t= 50)
    fig.show()

Order 1 (Forward Euler)¶

#####################################################
class ode_result:
    def __init__(self, y, t):
        self.y = y
        self.t = t

#####################################################
def forward_euler(tini, tend, nt, yini, fcn):

    dt = (tend-tini) / (nt-1)
    t = np.linspace(tini, tend, nt)

    yini = np.atleast_1d(yini)
    neq = yini.size

    y = np.zeros((neq, nt))
    y[:,0] = yini

    for it, tn  in enumerate(t[:-1]):
        yn = y[:,it]
        y[:,it+1] = yn + dt*np.array(fcn(tn, yn))

    return ode_result(y, t)

#####################################################
def stabiliy_disk_rk1():
    x = np.linspace(-3.5, 1, 500)
    y = np.linspace(-2.5, 2, 500)

    z = x + 1j*y[:, np.newaxis]
    disk = np.zeros_like(z, dtype = np.double)
    rk1 = z + 1
    mask = np.abs(rk1)<=1
    disk[mask] = np.abs(rk1[mask])
    return x, y, disk

uini = 2.
tini = 0.
tend = 1.5
k = 50.
plot_sol_and_error(uini, tini, tend, k, order=1)

¶

Order 2 (RK2)¶

#####################################################
def rk2(tini, tend, nt, yini, fcn):

    dt = (tend-tini) / (nt-1)
    t = np.linspace(tini, tend, nt)

    yini = np.atleast_1d(yini)
    neq = yini.size

    y = np.zeros((neq, nt))
    y[:,0] = yini

    for it, tn  in enumerate(t[:-1]):
        yn = y[:,it]
        k1 = fcn(tn, yn)
        k2 = fcn(tn + 0.5*dt, yn + dt*(0.5*k1))
        y[:,it+1] = yn + dt*k2

    return ode_result(y, t)

#####################################################
def stabiliy_disk_rk2():
    x = np.linspace(-3.5, 1.2, 500)
    y = np.linspace(-2.5, 2.2, 500)

    z = x + 1j*y[:, np.newaxis]
    disk = np.zeros_like(z, dtype = np.double)
    rk2 = z + z**2/2 + 1
    mask = np.abs(rk2)<=1
    disk[mask] = np.abs(rk2[mask])
    return x, y, disk

uini = 2.
tini = 0.
tend = 1.5
k = 50.
plot_sol_and_error(uini, tini, tend, k, order=2)

Ordre 3 (RK3)¶

#####################################################
def rk3(tini, tend, nt, yini, fcn):

    dt = (tend-tini) / (nt-1)
    t = np.linspace(tini, tend, nt)

    yini = np.atleast_1d(yini)
    neq = yini.size

    y = np.zeros((neq, nt))
    y[:,0] = yini

    for it, tn  in enumerate(t[:-1]):
        yn = y[:,it]
        k1 = fcn(tn, yn)
        k2 = fcn(tn + 0.5*dt, yn + dt*(0.5*k1))
        k3 = fcn(tn + dt, yn + dt*(-k1 + 2*k2))
        y[:,it+1] = yn + (dt/6)*(k1+4*k2+k3)

    return ode_result(y, t)

#####################################################
def stabiliy_disk_rk3():
    x = np.linspace(-4, 2, 500)
    y = np.linspace(-3, 3, 500)

    z = x + 1j*y[:, np.newaxis]
    disk = np.zeros_like(z, dtype = np.double)
    rk3 = z + z**2/2 + z**3/6 + 1
    mask = np.abs(rk3)<=1
    disk[mask] = np.abs(rk3[mask])
    return x, y, disk

uini = 2.
tini = 0.
tend = 1.5
k = 50.
plot_sol_and_error(uini, tini, tend, k, order=3)

Ordre 4 (RK4)¶

#####################################################
def rk4(tini, tend, nt, yini, fcn):

    dt = (tend-tini) / (nt-1)
    t = np.linspace(tini, tend, nt)

    yini = np.atleast_1d(yini)
    neq = yini.size

    y = np.zeros((neq, nt))
    y[:,0] = yini

    for it, tn  in enumerate(t[:-1]):
        yn = y[:,it]
        k1 = fcn(tn, yn)
        k2 = fcn(tn + 0.5*dt, yn + dt*(0.5*k1))
        k3 = fcn(tn + 0.5*dt, yn + dt*(0.5*k2))
        k4 = fcn(tn + dt, yn + dt*k3)
        y[:,it+1] = yn + (dt/6)*(k1+2*k2+2*k3+k4)

    return ode_result(y, t)

#####################################################
def stabiliy_disk_rk4():
    x = np.linspace(-4.5, 2.5, 500)
    y = np.linspace(-3.5, 3.5, 500)

    z = x + 1j*y[:, np.newaxis]
    disk = np.zeros_like(z, dtype = np.double)
    rk4 = z + z**2/2 + z**3/6 + z**4/24 + 1
    mask = np.abs(rk4)<=1
    disk[mask] = np.abs(rk4[mask])
    return x, y, disk

uini = 2.
tini = 0.
tend = 1.5
k = 50.
plot_sol_and_error(uini, tini, tend, k, order=4)