Notebook

Operator splitting methods¶

Operator splitting is a natural and old idea. When a PDE or system of PDEs contains different terms expressing different physics, it is natural to use different numerical methods for different physical processes. This can optimize and simplify the overall solution process. The idea was especially popularized in the context of the Navier-Stokes equations and reaction-diffusion PDEs. Common names for the technique are operator splitting, fractional step methods, and split-step methods. We shall stick to the former name. In the context of nonlinear differential equations, operator splitting can be used to isolate nonlinear terms and simplify the solution methods.

A related technique, often known as dimensional splitting or alternating direction implicit (ADI) methods, is to split the spatial dimensions and solve a 2D or 3D problem as two or three consecutive 1D problems, but this type of splitting is not to be further considered here.

Ordinary operator splitting for ODEs¶

Consider first an ODE where the right-hand side is split into two terms:

$$ \begin{equation} u' = f_0(u) + f_1(u)\thinspace . \label{_auto1} \tag{1} \end{equation} $$

In case $f_0$ and $f_1$ are linear functions of $u$, $f_0=au$ and $f_1=bu$, we have $u(t)=Ie^{(a+b)t}$, if $u(0)=I$. When going one time step of length $\Delta t$ from $t_n$ to $t_{n+1}$, we have

$$ u(t_{n+1}) = u(t_n)e^{(a+b)\Delta t}\thinspace . $$

This expression can be also be written as

$$ u(t_{n+1}) = u(t_n)e^{a\Delta t}e^{b\Delta t}, $$

$$ \begin{equation} u^\stepzero = u(t_n)e^{a\Delta t}, \label{nonlin:splitting:ODE:step1} \tag{2} \end{equation} $$

$$ \begin{equation} u(t_{n+1}) = u^\stepzero e^{b\Delta t} \label{nonlin:splitting:ODE:step2} \tag{3} \end{equation} $$

The first step (2) means solving $u'=f_0$ over a time interval $\Delta t$ with $u(t_n)$ as start value. The second step (3) means solving $u'=f_1$ over a time interval $\Delta t$ with the value at the end of the first step as start value. That is, we progress the solution in two steps and solve two ODEs $u'=f_0$ and $u'=f_1$. The order of the equations is not important. From the derivation above we see that solving $u'=f_1$ prior to $u'=f_0$ can equally well be done.

The technique is exact if the ODEs are linear. For nonlinear ODEs it is only an approximate method with error $\Delta t$. The technique can be extended to an arbitrary number of steps; i.e., we may split the PDE system into any number of subsystems. Examples will illuminate this principle.

Strang splitting for ODEs¶

The accuracy of the splitting method in the section Ordinary operator splitting for ODEs can be improved from $\Oof{\Delta t}$ to $\Oof{\Delta t^2}$ using so-called Strang splitting, where we take frac{1}{2} a step with the $f_0$ operator, a full step with the $f_1$ operator, and finally frac{1}{2} another step with the $f_0$ operator. During a time interval $\Delta t$ the algorithm can be written as follows.

$$ \begin{align*} \frac{du^{\stepzero}}{dt} &= f_0(u^{\stepzero}), \quad u^{\stepzero}(t_n)=u(t_n), \quad t\in [t_n,t_n+\frac{1}{2}\Delta t],\\ \frac{du^{\stephalf}}{dt} &= f_1(u^{\stephalf}), \quad u^{\stephalf}(t_n)=u^{\stepzero}(t_{n+\frac{1}{2}}), \quad t\in [t_n,t_n+\Delta t],\\ \frac{du^{\stepone}}{dt} &= f_0(u^{\stepone}), \quad u^{\stepone}(t_{n+\frac{1}{2}})=u^{\stephalf}(t_{n+1}), \quad t\in [t_n+\frac{1}{2}\Delta t, t_n+\Delta t]\thinspace . \end{align*} $$

The global solution is set as $u(t_{n+1}) = u^{\stepone}(t_{n+1})$.

There is no use in combining higher-order methods with ordinary splitting since the error due to splitting is $\Oof{\Delta t}$, but for Strang splitting it makes sense to use schemes of order $\Oof{\Delta t^2}$.

With the notation introduced for Strang splitting, we may express ordinary first-order splitting as

$$ \begin{align*} \frac{du^{\stepzero}}{dt} &= f_0(u^{\stepzero}),\quad u^{\stepzero}(t_n)=u(t_n), \quad t\in [t_n,t_n+\Delta t],\\ \frac{du^{\stepone}}{dt} &= f_1(u^{\stepone}),\quad u^{\stepone}(t_n)=u^{\stepzero}(t_{n+1}), \quad t\in [t_n,t_n+\Delta t], \end{align*} $$

with global solution set as $u(t_{n+1}) = u^{\stepone}(t_{n+1})$.

Example: Logistic growth¶

Let us split the (scaled) logistic equation

$$ u'=u(1-u),\quad u(0)=0.1, $$

with solution $u=(9e^{-t}+1)^{-1}$, into

$$ u'=u - u^2 = f_0(u) + f_1(u), \quad f_0(u)=u,\quad f_1(u)=-u^2\thinspace . $$

We solve $u'=f_0(u)$ and $u'=f_1(u)$ by a Forward Euler step. In addition, we add a method where we solve $u'=f_0(u)$ analytically, since the equation is actually $u'=u$ with solution $e^t$. The software that accompanies the following methods is the file split_logistic.py.

Splitting techniques¶

Ordinary splitting takes a Forward Euler step for each of the ODEs according to

$$ \begin{equation} \frac{u^{\stepzero,n+1} - u^{\stepzero,n}}{\Delta t} = f_0(u^{\stepzero,n}),\quad u^{\stepzero,n} = u(t_n),\quad t\in [t_n,t_n+\Delta t], \label{_auto2} \tag{4} \end{equation} $$

$$ \begin{equation} \frac{u^{\stepone,n+1} - u^{\stepone, n}}{\Delta t} = f_1(u^{\stepone,n}),\quad u^{\stepone,n} = u^{\stepzero,n+1},\quad t\in [t_n,t_n+\Delta t], \label{_auto3} \tag{5} \end{equation} $$

with $u(t_{n+1}) = u^{\stepone,n+1}$.

Strang splitting takes the form

$$ \begin{equation} \frac{u^{\stepzero,n+\frac{1}{2}} - u^{\stepzero,n}}{\frac{1}{2}\Delta t} = f_0(u^{\stepzero,n}),\quad u^{\stepzero,n} = u(t_n),\quad t\in [t_n,t_n+\frac{1}{2}\Delta t], \label{_auto4} \tag{6} \end{equation} $$

$$ \begin{equation} \frac{u^{\stephalf,n+1}-u^{\stephalf,n}}{\Delta t} = f_1(u^{\stephalf,n}),\quad u^{\stephalf,n} = u^{\stepzero, n+\frac{1}{2}},\quad t\in [t_n,t_n+\Delta t], \label{_auto5} \tag{7} \end{equation} $$

$$ \begin{equation} \frac{u^{\stepone, n+1} - u^{\stepone, n+\frac{1}{2}}}{\frac{1}{2}\Delta t} = f_0(u^{\stepone,n+\frac{1}{2}}),\quad u^{\stepone,n+\frac{1}{2}} = u^{\stephalf,n+1},\quad t\in [t_n+\frac{1}{2}\Delta t, t_n+\Delta t]\thinspace . \label{_auto6} \tag{8} \end{equation} $$

Verbose implementation¶

The following function computes four solutions arising from the Forward Euler method, ordinary splitting, Strang splitting, as well as Strang splitting with exact treatment of $u'=f_0(u)$:

In [1]:

import numpy as np

def solver(dt, T, f, f_0, f_1):
    """
    Solve u'=f by the Forward Euler method and by ordinary and
    Strang splitting: f(u) = f_0(u) + f_1(u).
    """
    Nt = int(round(T/float(dt)))
    t = np.linspace(0, Nt*dt, Nt+1)
    u_FE = np.zeros(len(t))
    u_split1 = np.zeros(len(t))  # 1st-order splitting
    u_split2 = np.zeros(len(t))  # 2nd-order splitting
    u_split3 = np.zeros(len(t))  # 2nd-order splitting w/exact f_0

    # Set initial values
    u_FE[0] = 0.1
    u_split1[0] = 0.1
    u_split2[0] = 0.1
    u_split3[0] = 0.1

    for n in range(len(t)-1):
        # Forward Euler method
        u_FE[n+1] = u_FE[n] + dt*f(u_FE[n])

        # --- Ordinary splitting ---
        # First step
        u_s_n = u_split1[n]
        u_s = u_s_n + dt*f_0(u_s_n)
        # Second step
        u_ss_n = u_s
        u_ss = u_ss_n + dt*f_1(u_ss_n)
        u_split1[n+1] = u_ss

        # --- Strang splitting ---
        # First step
        u_s_n = u_split2[n]
        u_s = u_s_n + dt/2.*f_0(u_s_n)
        # Second step
        u_sss_n = u_s
        u_sss = u_sss_n + dt*f_1(u_sss_n)
        # Third step
        u_ss_n = u_sss
        u_ss = u_ss_n + dt/2.*f_0(u_ss_n)
        u_split2[n+1] = u_ss

        # --- Strang splitting using exact integrator for u'=f_0 ---
        # First step
        u_s_n = u_split3[n]
        u_s = u_s_n*np.exp(dt/2.)  # exact
        # Second step
        u_sss_n = u_s
        u_sss = u_sss_n + dt*f_1(u_sss_n)
        # Third step
        u_ss_n = u_sss
        u_ss = u_ss_n*np.exp(dt/2.)  # exact
        u_split3[n+1] = u_ss

    return u_FE, u_split1, u_split2, u_split3, t

Compact implementation¶

We have used quite many lines for the steps in the splitting methods. Many will prefer to condense the code a bit, as done here:

In [2]:

# Ordinary splitting
u_s = u_split1[n] + dt*f_0(u_split1[n])
u_split1[n+1] = u_s + dt*f_1(u_s)

# Strang splitting
u_s = u_split2[n] + dt/2.*f_0(u_split2[n])
u_sss = u_s + dt*f_1(u_s)
u_split2[n+1] = u_sss + dt/2.*f_0(u_sss)

# Strang splitting using exact integrator for u'=f_0
u_s = u_split3[n]*np.exp(dt/2.)  # exact
u_ss = u_s + dt*f_1(u_s)
u_split3[n+1] = u_ss*np.exp(dt/2.)

Results¶

Figure shows that the impact of splitting is significant. Interestingly, however, the Forward Euler method applied to the entire problem directly is much more accurate than any of the splitting schemes. We also see that Strang splitting is definitely more accurate than ordinary splitting and that it helps a bit to use an exact solution of $u'=f_0(u)$. With a large time step ($\Delta t = 0.2$, left plot in Figure), the asymptotic values are off by 20-30%. A more reasonable time step ($\Delta t = 0.05$, right plot in Figure) gives better results, but still the asymptotic values are up to 10% wrong.

As technique for solving nonlinear ODEs, we realize that the present case study is not particularly promising, as the Forward Euler method both linearizes the original problem and provides a solution that is much more accurate than any of the splitting techniques. In complicated multi-physics settings, on the other hand, splitting may be the only feasible way to go, and sometimes you really need to apply different numerics to different parts of a PDE problem. But in very simple problems, like the logistic ODE, splitting is just an inferior technique. Still, the logistic ODE is ideal for introducing all the mathematical details and for investigating the behavior.

Effect of ordinary and Strang splitting for the logistic equation.

Reaction-diffusion equation¶

Consider a diffusion equation coupled to chemical reactions modeled by a nonlinear term $f(u)$:

$$ \frac{\partial u}{\partial t} = \dfc\nabla^2u + f(u)\thinspace . $$

This is a physical process composed of two individual processes: $u$ is the concentration of a substance that is locally generated by a chemical reaction $f(u)$, while $u$ is spreading in space because of diffusion. There are obviously two time scales: one for the chemical reaction and one for diffusion. Typically, fast chemical reactions require much finer time stepping than slower diffusion processes. mathcal{I}_t could therefore be advantageous to split the two physical effects in separate models and use different numerical methods for the two.

A natural spitting in the present case is

$$ \begin{equation} \frac{\partial u^{\stepzero}}{\partial t} = \dfc\nabla^2 u^{\stepzero}, \label{nonlin:splitting:RD:eq:diffu} \tag{9} \end{equation} $$

$$ \begin{equation} \frac{\partial u^{\stepone}}{\partial t} = f(u^{\stepone}) \label{nonlin:splitting:RD:eq:decay} \tag{10}\thinspace . \end{equation} $$

Looking at these familiar problems, we may apply a $\theta$ rule (implicit) scheme for (9) over one time step and avoid dealing with nonlinearities by applying an explicit scheme for (10) over the same time step.

Suppose we have some solution $u$ at time level $t_n$. For flexibility, we define a $\theta$ method for the diffusion part (9) by

$$ [D_t u^\stepzero = \dfc (D_xD_x u^\stepzero + D_y D_y u^\stepzero)]^{n+\theta}\thinspace . $$

We use $u^{n}$ as initial condition for $u^\stepzero$.

The reaction part, which is defined at each mesh point (without coupling values in different mesh points), can employ any scheme for an ODE. Here we use an Adams-Bashforth method of order 2. Recall that the overall accuracy of the splitting method is maximum $\Oof{\Delta t^2}$ for Strang splitting, otherwise it is just $\Oof{\Delta t}$. Higher-order methods for ODEs will therefore be a waste of work. The 2nd-order Adams-Bashforth method reads

$$ \begin{equation} u^{\stepone,n+1}_{i,j} = u^{\stepone,n}_{i,j} + \frac{1}{2}\Delta t\left( 3f(u^{\stepone, n}_{i,j}, t_n) - f(u^{\stepone, n-1}_{i,j}, t_{n-1}) \right) \thinspace . \label{_auto7} \tag{11} \end{equation} $$

We can use a Forward Euler step to start the method, i.e, compute $u^{\stepone,1}_{i,j}$.

The algorithm goes like this:

Solve the diffusion problem for one time step as usual.
Solve the reaction ODEs at each mesh point in $[t_n,t_n+\Delta t]$, using the diffusion solution in 1. as initial condition. The solution of the ODEs constitutes the solution of the original problem at the end of each time step.

We may use a much smaller time step when solving the reaction part, adapted to the dynamics of the problem $u'=f(u)$. This gives great flexibility in splitting methods.

[hpl 1: Can you make an example? I have started in src/split_diffu_react.py.]

Example: Reaction-Diffusion with linear reaction term¶

The methods above may be explored in detail through a specific computational example in which we compute the convergence rates associated with four different solution approaches for the reaction-diffusion equation with a linear reaction term, i.e. $f(u)=-bu$. The methods comprise solving without splitting (just straight Forward Euler), ordinary splitting, first order Strang splitting, and second order Strang splitting. In all four methods, a standard centered difference approximation is used for the spatial second derivative. The methods share the error model $E = C h^r$, while differing in the step $h$ (being either $\Delta x^2$ or $\Delta x$) and the convergence rate $r$ (being either 1 or 2).

All code commented below is found in the file split_diffu_react.py. When executed, a function convergence_rates is called, from which all convergence rate computations are handled:

In [3]:

def convergence_rates(scheme='diffusion'):
    
    F = 0.5     # Upper limit for FE (stability). For CN, this 
                # limit does not apply, but for simplicity, we
                # choose F = 0.5 as the initial F value.
    T = 1.2
    a = 3.5
    b = 1
    L = 1.5
    k = np.pi/L
    
    def exact(x, t):
        '''exact sol. to: du/dt = a*d^2u/dx^2 - b*u'''
        return np.exp(-(a*k**2 + b)*t) * np.sin(k*x)
        
    def f(u, t):
        return -b*u  

    def I(x):
        return exact(x, 0)
        
    global error    # error computed in the user action function
    error = 0  

    # Convergence study
    def action(u, x, t, n):
        global error            
        if n == 1:      # New simulation, - reset error
            error = 0
        else:
            error = max(error, np.abs(u - exact(x, t[n])).max())        

    E = []
    h = []
    Nx_values = [10, 20, 40, 80]   # i.e., dx halved each time
    for Nx in Nx_values:        
        dx = L/Nx           
        if scheme == 'Strang_splitting_2ndOrder':
            print 'Strang splitting with 2nd order schemes...'
            # In this case, E = C*h**r (with r = 2) and since 
            # h = dx = K*dt, the ratio dt/dx must be constant.
            # To fulfill this demand, we must let F change
            # when dx changes. From F = a*dt/dx**2, it follows
            # that halving dx AND doubling F assures dt/dx const.
            # Initially, we simply choose F = 0.5.

            dt = F/a*dx**2
            #print 'dt/dx:', dt/dx            
            Nt = int(round(T/float(dt)))
            t = np.linspace(0, Nt*dt, Nt+1)   # global time            
            Strang_splitting_2ndOrder(I=I, a=a, b=b, f=f, L=L, dt=dt, 
                                      dt_Rfactor=1, F=F, t=t, T=T,
                                      user_action=action)    
            h.append(dx)
            # prepare for next iteration (make F match dx/2)
            F = F*2       # assures dt/dx const. when dx = dx/2                                      
        else:   
            # In these cases, E = C*h**r (with r = 1) and since 
            # h = dx**2 = K*dt, the ratio dt/dx**2 must be constant.
            # This is fulfilled by choosing F = 0.5 (for FE stability)
            # and make sure that F, dx and dt comply to F = a*dt/dx**2.            
            dt = F/a*dx**2
            Nt = int(round(T/float(dt)))
            t = np.linspace(0, Nt*dt, Nt+1)   # global time
            if scheme == 'diffusion':
                print 'FE on whole eqn...'
                diffusion_theta(I, a, f, L, dt, F, t, T,
                                step_no=0, theta=0,
                                u_L=0, u_R=0, user_action=action)                         
                h.append(dx**2)
            elif scheme == 'ordinary_splitting':
                print 'Ordinary splitting...'
                ordinary_splitting(I=I, a=a, b=b, f=f, L=L, dt=dt, 
                                   dt_Rfactor=1, F=F, t=t, T=T,
                                   user_action=action)        
                h.append(dx**2)
            elif scheme == 'Strang_splitting_1stOrder':
                print 'Strang splitting with 1st order schemes...'
                Strang_splitting_1stOrder(I=I, a=a, b=b, f=f, L=L, dt=dt, 
                                          dt_Rfactor=1, F=F, t=t, T=T,
                                          user_action=action)       
                h.append(dx**2)
            else:
                print 'Unknown scheme requested!'
                sys.exit(0)
        
            #print 'dt/dx**2:', dt/dx**2            
        
        E.append(error)
        Nx *= 2         # Nx doubled gives dx/2

    print 'E:', E
    print 'h:', h
        
    # Convergence rates 
    r = [np.log(E[i]/E[i-1])/np.log(h[i]/h[i-1])
         for i in range(1,len(Nx_values))]
    print 'Computed rates:', r


if __name__ == '__main__':
    
    schemes = ['diffusion', 
               'ordinary_splitting', 
               'Strang_splitting_1stOrder',
               'Strang_splitting_2ndOrder']
               
    for scheme in schemes:
        convergence_rates(scheme=scheme)

Now, with respect to the error ($E = C h^r$), the Forward Euler scheme, the ordinary splitting scheme and first order Strang splitting scheme are all first order ($r = 1$), with a step $h = \Delta x^2 = K^{-1}\Delta t$, where $K$ is some constant. This implies that the ratio $\frac{\Delta t}{\Delta x^2}$ must be held constant during convergence rate calculations. Furthermore, the Fourier number $F = \frac{\alpha \Delta t}{\Delta x^2}$ is upwards limited to $F = 0.5$, being the stability limit with explicit schemes. Thus, in these cases, we use the fixed value of $F$ and a given (but changing) spatial resolution $\Delta x$ to compute the corresponding value of $\Delta t$ according to the expression for $F$. This assures that $\frac{\Delta t}{\Delta x^2}$ is kept constant. The loop in convergence_rates runs over a chosen set of grid points (Nx_values) which gives a doubling of spatial resolution with each iteration ($\Delta x$ is halved).

For the second order Strang splitting scheme, we have $r = 2$ and a step $h = \Delta x = K^{-1}\Delta t$, where $K$ again is some constant. In this case, it is thus the ratio $\frac{\Delta t}{\Delta x}$ that must be held constant during the convergence rate calculations. From the expression for $F$, it is clear then that $F$ must change with each halving of $\Delta x$. In fact, if $F$ is doubled each time $\Delta x$ is halved, the ratio $\frac{\Delta t}{\Delta x}$ will be constant (this follows, e.g., from the expression for $F$). This is utilized in our code.

A solver diffusion_theta is used in each of the four solution approaches:

In [4]:

def diffusion_theta(I, a, f, L, dt, F, t, T, step_no, theta=0.5, 
                    u_L=0, u_R=0, user_action=None):
    """
    Full solver for the model problem using the theta-rule
    difference approximation in time (no restriction on F,
    i.e., the time step when theta >= 0.5).     Vectorized 
    implementation and sparse (tridiagonal) coefficient matrix.
    Note that t always covers the whole global time interval, whether
    splitting is the case or not. T, on the other hand, is 
    the end of the global time interval if there is no split,
    but if splitting, we use T=dt. When splitting, step_no 
    keeps track of the time step number (for lookup in t).    
    """

    Nt = int(round(T/float(dt)))
    dx = np.sqrt(a*dt/F)   
    Nx = int(round(L/dx))
    x = np.linspace(0, L, Nx+1)       # Mesh points in space
    # Make sure dx and dt are compatible with x and t
    dx = x[1] - x[0]
    dt = t[1] - t[0]

    u   = np.zeros(Nx+1)   # solution array at t[n+1]
    u_1 = np.zeros(Nx+1)   # solution at t[n]

    # Representation of sparse matrix and right-hand side
    diagonal = np.zeros(Nx+1)
    lower    = np.zeros(Nx)
    upper    = np.zeros(Nx)
    b        = np.zeros(Nx+1)

    # Precompute sparse matrix (scipy format)
    Fl = F*theta
    Fr = F*(1-theta)
    diagonal[:] = 1 + 2*Fl
    lower[:] = -Fl  #1
    upper[:] = -Fl  #1
    # Insert boundary conditions
    diagonal[0] = 1
    upper[0] = 0
    diagonal[Nx] = 1
    lower[-1] = 0

    diags = [0, -1, 1]
    A = scipy.sparse.diags(
        diagonals=[diagonal, lower, upper],
        offsets=[0, -1, 1], shape=(Nx+1, Nx+1),
        format='csr')
    #print A.todense()

    # Allow f to be None or 0
    if f is None or f == 0:
        f = lambda x, t: np.zeros((x.size)) \
            if isinstance(x, np.ndarray) else 0

    # Set initial condition  
    if isinstance(I, np.ndarray):   # I is an array
        u_1 = np.copy(I)
    else:                           # I is a function
        for i in range(0, Nx+1):
            u_1[i] = I(x[i])

    if user_action is not None:
        user_action(u_1, x, t, step_no+0)        

    # Time loop
    for n in range(0, Nt):
        b[1:-1] = u_1[1:-1] + \
                  Fr*(u_1[:-2] - 2*u_1[1:-1] + u_1[2:]) + \
                  dt*theta*f(u_1[1:-1], t[step_no+n+1]) + \
                  dt*(1-theta)*f(u_1[1:-1], t[step_no+n])
        b[0] = u_L; b[-1] = u_R  # boundary conditions
        u[:] = scipy.sparse.linalg.spsolve(A, b)

        if user_action is not None:
            user_action(u, x, t, step_no+(n+1))

        # Update u_1 before next step
        u_1, u = u, u_1

    # u is now contained in u_1 (swapping)
    return u_1

For the no splitting approach with Forward Euler in time, this solver handles both the diffusion and the reaction term. When splitting, diffusion_theta takes care of the diffusion term only, while the reaction term is handled either by a Forward Euler scheme in reaction_FE, or by a second order Adams-Bashforth scheme from Odespy. The reaction_FE function covers one complete time step dt during ordinary splitting, while Strang splitting (both first and second order) applies it with dt/2 twice during each time step dt. Since the reaction term typically represents a much faster process than the diffusion term, a further refinement of the time step is made possible in reaction_FE. mathcal{I}_t was implemented as

In [5]:

def reaction_FE(I, f, L, Nx, dt, dt_Rfactor, t, step_no, 
                user_action=None):
    """Reaction solver, Forward Euler method.
    Note the at t covers the whole global time interval. 
    dt is either one complete,or one frac{1}{2}, of the step in the 
    diffusion part, i.e. there is a local time interval 
    [0, dt] or [0, dt/2] that the reaction_FE
    deals with each time it is called. step_no keeps
    track of the (global) time step number (required 
    for lookup in t).    
    """

    u = np.copy(I)      
    dt_local = dt/float(dt_Rfactor)    
    Nt_local = int(round(dt/float(dt_local)))    
    x = np.linspace(0, L, Nx+1)  
          
    for n in range(Nt_local):
        time = t[step_no] + n*dt_local
        u[1:Nx] = u[1:Nx] + dt_local*f(u[1:Nx], time)  
        
    # BC already inserted in diffusion step, i.e. no action here            
    return u

With the ordinary splitting approach, each time step dt is covered twice. First computing the impact of the reaction term, then the contribution from the diffusion term:

In [6]:

def ordinary_splitting(I, a, b, f, L, dt, 
                       dt_Rfactor, F, t, T, 
                       user_action=None):
    '''1st order scheme, i.e. Forward Euler is enough for both
    the diffusion and the reaction part. The time step dt is 
    given for the diffusion step, while the time step for the 
    reaction part is found as dt/dt_Rfactor, where dt_Rfactor >= 1.
    '''
    Nt = int(round(T/float(dt)))
    dx = np.sqrt(a*dt/F)
    Nx = int(round(L/dx))
    x = np.linspace(0, L, Nx+1)       # Mesh points in space
    u = np.zeros(Nx+1)
        
    # Set initial condition u(x,0) = I(x)
    for i in range(0, Nx+1):
        u[i] = I(x[i])

    # In the following loop, each time step is "covered twice", 
    # first for reaction, then for diffusion
    for n in range(0, Nt):    
        # Reaction step (potentially many smaller steps within dt)
        u_s = reaction_FE(I=u, f=f, L=L, Nx=Nx, 
                        dt=dt, dt_Rfactor=dt_Rfactor, 
                        t=t, step_no=n, 
                        user_action=None)        
        
        u = diffusion_theta(I=u_s, a=a, f=0, L=L, dt=dt, F=F,
                              t=t, T=dt, step_no=n, theta=0,
                              u_L=0, u_R=0, user_action=None)                         
                           
        if user_action is not None:
            user_action(u, x, t, n+1)

    return

For the two Strang splitting approaches, each time step dt is handled by first computing the reaction step for (the first) dt/2, followed by a diffusion step dt, before the reaction step is treated once again for (the remaining) dt/2. Since first order Strang splitting is no better than first order accurate, both the reaction and diffusion steps are computed explicitly. The solver was implemented as

In [7]:

def Strang_splitting_1stOrder(I, a, b, f, L, dt, dt_Rfactor, 
                              F, t, T, user_action=None):
    '''Strang splitting while still using FE for the reaction
    step and for the diffusion step. Gives 1st order scheme.
    The time step dt is given for the diffusion step, while 
    the time step for the reaction part is found as 
    0.5*dt/dt_Rfactor, where dt_Rfactor >= 1. Introduce an 
    extra time mesh t2 for the reaction part, since it steps dt/2.
    '''
    Nt = int(round(T/float(dt)))
    t2 = np.linspace(0, Nt*dt, (Nt+1)+Nt)   # Mesh points in diff    
    dx = np.sqrt(a*dt/F)
    Nx = int(round(L/dx))
    x = np.linspace(0, L, Nx+1)     
    u = np.zeros(Nx+1)
        
    # Set initial condition u(x,0) = I(x)
    for i in range(0, Nx+1):
        u[i] = I(x[i])

    for n in range(0, Nt):    
        # Reaction step (1/2 dt: from t_n to t_n+1/2)
        # (potentially many smaller steps within dt/2)
        u_s = reaction_FE(I=u, f=f, L=L, Nx=Nx, 
                        dt=dt/2.0, dt_Rfactor=dt_Rfactor, 
                        t=t2, step_no=2*n, 
                        user_action=None)
        # Diffusion step (1 dt: from t_n to t_n+1)
        u_sss = diffusion_theta(I=u_s, a=a, f=0, L=L, dt=dt, F=F,
                              t=t, T=dt, step_no=n, theta=0,
                              u_L=0, u_R=0, user_action=None)                         
        # Reaction step (1/2 dt: from t_n+1/2 to t_n+1)
        # (potentially many smaller steps within dt/2)
        u = reaction_FE(I=u_sss, f=f, L=L, Nx=Nx, 
                        dt=dt/2.0, dt_Rfactor=dt_Rfactor, 
                        t=t2, step_no=2*n+1, 
                        user_action=None)
        
        if user_action is not None:
            user_action(u, x, t, n+1)

    return

The second order version of the Strang splitting approach utilizes a second order Adams-Bashforth solver for the reaction part and a Crank-Nicolson scheme for the diffusion part. The solver has the same structure as the one for first order Strang splitting and was implemented as

In [8]:

def Strang_splitting_2ndOrder(I, a, b, f, L, dt, dt_Rfactor, 
                              F, t, T, user_action=None):
    '''Strang splitting using Crank-Nicolson for the diffusion
    step (theta-rule) and Adams-Bashforth 2 for the reaction step. 
    Gives 2nd order scheme. Introduce an extra time mesh t2 for 
    the reaction part, since it steps dt/2.
    '''
    import odespy
    Nt = int(round(T/float(dt)))
    t2 = np.linspace(0, Nt*dt, (Nt+1)+Nt)   # Mesh points in diff    
    dx = np.sqrt(a*dt/F)
    Nx = int(round(L/dx))
    x = np.linspace(0, L, Nx+1)     
    u = np.zeros(Nx+1)
        
    # Set initial condition u(x,0) = I(x)
    for i in range(0, Nx+1):
        u[i] = I(x[i])

    reaction_solver = odespy.AdamsBashforth2(f)    

    for n in range(0, Nt):    
        # Reaction step (1/2 dt: from t_n to t_n+1/2)
        # (potentially many smaller steps within dt/2)
        reaction_solver.set_initial_condition(u)
        t_points = np.linspace(0, dt/2.0, dt_Rfactor+1)
        u_AB2, t_ = reaction_solver.solve(t_points) # t_ not needed
        u_s = u_AB2[-1,:]  # pick sol at last point in time

        # Diffusion step (1 dt: from t_n to t_n+1)
        u_sss = diffusion_theta(I=u_s, a=a, f=0, L=L, dt=dt, F=F,
                              t=t, T=dt, step_no=n, theta=0.5,
                              u_L=0, u_R=0, user_action=None)                         
        # Reaction step (1/2 dt: from t_n+1/2 to t_n+1)
        # (potentially many smaller steps within dt/2)
        reaction_solver.set_initial_condition(u_sss)
        t_points = np.linspace(0, dt/2.0, dt_Rfactor+1)
        u_AB2, t_ = reaction_solver.solve(t_points) # t_ not needed
        u = u_AB2[-1,:]  # pick sol at last point in time
               
        if user_action is not None:
            user_action(u, x, t, n+1)

    return

When executing split_diffu_react.py, we find that the estimated convergence rates are as expected. The second order Strang splitting gives the least error (about $4e^{-5}$) and has second order convergence ($r = 2$), while the remaining three approaches have first order convergence ($r = 1$).

Analysis of the splitting method¶

Let us address a linear PDE problem for which we can develop analytical solutions of the discrete equations, with and without splitting, and discuss these. Choosing $f(u)=-\beta u$ for a constant $\beta$ gives a linear problem. We use the Forward Euler method for both the PDE and ODE problems.

We seek a 1D Fourier wave component solution of the problem, assuming homogeneous Dirichlet conditions at $x=0$ and $x=L$:

$$ u = e^{-\dfc k^2 t - \beta t}\sin kx,\quad k = \frac{\pi}{L}\thinspace . $$

This component fits the 1D PDE problem ($f=0$). On complex form we can write

$$ u = e^{-\dfc k^2 t - \beta t + ikx}, $$

where $i=\sqrt{-1}$ and the imaginary part is taken as the physical solution.

We refer to the section "Analysis of schemes for the diffusion equation": "" in [Langtangen_Linge] and to the book [Langtangen_decay] for a discussion of exact numerical solutions to diffusion and decay problems, respectively. The key idea is to search for solutions $A^ne^{ikx}$ and determine $A$. For the diffusion problem solved by a Forward Euler method one has

$$ A = 1 - 4F\sin^p, $$

where $F=\dfc\Delta t/\Delta x^2$ is the mesh Fourier number and $p=k\Delta x/2$ is a dimensionless number reflecting the spatial resolution (number of points per wave length in space). For the decay problem $u'=-\beta u$, we have $A=1 - q$, where $q$ is a dimensionless parameter for the resolution in the decay problem: $q = \beta\Delta t$.

The original model problem can also be discretized by a Forward Euler scheme,

$$ [D^+_t u = \dfc D_xD_x u - \beta u]^n_i\thinspace . $$

Assuming $A^ne^{ikx}$ we find that

$$ u^n_i = (1 - 4F\sin^p -q)^n\sin kx\thinspace . $$

We are particularly interested in what happens at one time step. That is,

$$ u^{n}_{i} = (1-4F\sin^2 p)u^{n-1}_i\thinspace . $$

In the two stage algorithm, we first compute the diffusion step

$$ u^{\stepzero,n+1}_i = (1 - 4F\sin^2 p)u^{n-1}_i\thinspace . $$

Then we use this as input to the decay algorithm and arrive at

$$ u^{\stepone,n+1} = (1-q)u^{\stepzero,n+1} = (1-q)(1-4F\sin^2 p) u^{n-1}_i\thinspace . $$

The splitting approximation over one step is therefore

$$ E = 1 - 4F\sin^p -q - (1-q)(1-4F\sin^2 p) = -q(2 - F\sin^2 p))\thinspace . $$