.. !split .. _nonlin:splitting: 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. .. _nonlin:splitting:ODE: Ordinary operator splitting for ODEs ------------------------------------ .. index:: operator splitting .. index:: splitting ODEs .. index:: fractional step methods Consider first an ODE where the right-hand side is split into two terms: .. _Eq:_auto245: .. math:: \tag{590} u' = f_0(u) + f_1(u){\thinspace .} In case :math:`f_0` and :math:`f_1` are linear functions of :math:`u`, :math:`f_0=au` and :math:`f_1=bu`, we have :math:`u(t)=Ie^{(a+b)t}`, if :math:`u(0)=I`. When going one time step of length :math:`\Delta t` from :math:`t_n` to :math:`t_{n+1}`, we have .. math:: u(t_{n+1}) = u(t_n)e^{(a+b)\Delta t}{\thinspace .} This expression can be also be written as .. math:: u(t_{n+1}) = u(t_n)e^{a\Delta t}e^{b\Delta t}, or .. _Eq:nonlin:splitting:ODE:step1: .. math:: \tag{591} u^{*} = u(t_n)e^{a\Delta t}, .. _Eq:nonlin:splitting:ODE:step2: .. math:: \tag{592} u(t_{n+1}) = u^{*} e^{b\Delta t} The first step :ref:`(591) ` means solving :math:`u'=f_0` over a time interval :math:`\Delta t` with :math:`u(t_n)` as start value. The second step :ref:`(592) ` means solving :math:`u'=f_1` over a time interval :math:`\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 :math:`u'=f_0` and :math:`u'=f_1`. The order of the equations is not important. From the derivation above we see that solving :math:`u'=f_1` prior to :math:`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 :math:`\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. .. _nonlin:splitting:ODE_Strang: Strang splitting for ODEs ------------------------- .. index:: Strang splitting The accuracy of the splitting method in the section :ref:`nonlin:splitting:ODE` can be improved from :math:`{\mathcal{O}(\Delta t)}` to :math:`{\mathcal{O}(\Delta t^2)}` using so-called *Strang splitting*, where we take half a step with the :math:`f_0` operator, a full step with the :math:`f_1` operator, and finally half another step with the :math:`f_0` operator. During a time interval :math:`\Delta t` the algorithm can be written as follows. .. math:: \begin{align*} \frac{du^{{*}}}{dt} &= f_0(u^{{*}}), \quad u^{{*}}(t_n)=u(t_n), \quad t\in [t_n,t_n+\frac{1}{2}\Delta t],\\ \frac{du^{{***}}}{dt} &= f_1(u^{{***}}), \quad u^{{***}}(t_n)=u^{{*}}(t_{n+\frac{1}{2}}), \quad t\in [t_n,t_n+\Delta t],\\ \frac{du^{{**}}}{dt} &= f_0(u^{{**}}), \quad u^{{**}}(t_n+\frac{1}{2})=u^{{***}}(t_{n+\frac{1}{2}}), \quad t\in [t_n+\frac{1}{2}\Delta t, t_n+\Delta t]{\thinspace .} \end{align*} The global solution is set as :math:`u(t_{n+1}) = u^{{**}}(t_{n+1})`. There is no use in combining higher-order methods with ordinary splitting since the error due to splitting is :math:`\Oof{\Delta t}`, but for Strang splitting it makes sense to use schemes of order :math:`{\mathcal{O}(\Delta t^2)}`. With the notation introduced for Strang splitting, we may express ordinary first-order splitting as .. math:: \begin{align*} \frac{du^{{*}}}{dt} &= f_0(u^{{*}}),\quad u^{{*}}(t_n)=u(t_n), \quad t\in [t_n,t_n+\Delta t],\\ \frac{du^{{**}}}{dt} &= f_1(u^{{**}}),\quad u^{{**}}(t_n)=u^{{*}}(t_{n+1}), \quad t\in [t_n,t_n+\Delta t], \end{align*} with global solution set as :math:`u(t_{n+1}) = u^{{**}}(t_{n+1})`. .. _nonlin:splitting:logistic: Example: Logistic growth ------------------------ Let us split the (scaled) logistic equation .. math:: u'=u(1-u),\quad u(0)=0.1, with solution :math:`u=(9e^{-t}+1)^{-1}`, into .. math:: u'=u - u^2 = f_0(u) + f_1(u), \quad f_0(u)=u,\quad f_1(u)=-u^2{\thinspace .} We solve :math:`u'=f_0(u)` and :math:`u'=f_1(u)` by a Forward Euler step. In addition, we add a method where we solve :math:`u'=f_0(u)` analytically, since the equation is actually :math:`u'=u` with solution :math:`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 .. _Eq:_auto246: .. math:: \tag{593} \frac{u^{{*},n+1} - u^{{*},n}}{\Delta t} = f_0(u^{{*},n}),\quad u^{{*},n} = u(t_n),\quad t\in [t_n,t_n+\Delta t], .. _Eq:_auto247: .. math:: \tag{594} \frac{u^{{**},n+1} - u^{{**}, n}}{\Delta t} = f_1(u^{{**},n}),\quad u^{{**},n} = u^{{*},n+1},\quad t\in [t_n,t_n+\Delta t], with :math:`u(t_{n+1}) = u^{{**},n+1}`. Strang splitting takes the form .. _Eq:_auto248: .. math:: \tag{595} \frac{u^{{*},n+\frac{1}{2}} - u^{{*},n}}{\frac{1}{2}\Delta t} = f_0(u^{{*},n}),\quad u^{{*},n} = u(t_n),\quad t\in [t_n,t_n+\frac{1}{2}\Delta t], .. _Eq:_auto249: .. math:: \tag{596} \frac{u^{{***},n+1}-u^{{***},n}}{\Delta t} = f_1(u^{{***},n}),\quad u^{{***},n} = u^{{*}, n+\frac{1}{2}},\quad t\in [t_n,t_n+\Delta t], .. _Eq:_auto250: .. math:: \tag{597} \frac{u^{{**}, n+1} - u^{{**}, n+\frac{1}{2}}}{\frac{1}{2}\Delta t} = f_0(u^{{**},n+\frac{1}{2}}),\quad u^{{**},n+\frac{1}{2}} = u^{{***},n+1},\quad t\in [t_n+\frac{1}{2}\Delta t, t_n+\Delta t]{\thinspace .} 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 :math:`u'=f_0(u)`: .. code-block:: python 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: .. code-block:: python # 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 :ref:`nonlin:splitting:ODE:logistic:fig` 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 :math:`u'=f_0(u)`. With a large time step (:math:`\Delta t = 0.2`, left plot in Figure :ref:`nonlin:splitting:ODE:logistic:fig`), the asymptotic values are off by 20-30%. A more reasonable time step (:math:`\Delta t = 0.05`, right plot in Figure :ref:`nonlin:splitting:ODE:logistic:fig`) 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. .. _nonlin:splitting:ODE:logistic:fig: .. figure:: split_logistic.png :width: 800 *Effect of ordinary and Strang splitting for the logistic equation* .. Cannot analyze logistic method because we have no exact solution .. of the discrete equations. But in the next linear problem we can analyze. .. _nonlin:splitting:RD: Reaction-diffusion equation --------------------------- Consider a diffusion equation coupled to chemical reactions modeled by a nonlinear term :math:`f(u)`: .. math:: \frac{\partial u}{\partial t} = {\alpha}\nabla^2u + f(u){\thinspace .} This is a physical process composed of two individual processes: :math:`u` is the concentration of a substance that is locally generated by a chemical reaction :math:`f(u)`, while :math:`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. It 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 .. _Eq:nonlin:splitting:RD:eq:diffu: .. math:: \tag{598} \frac{\partial u^{{*}}}{\partial t} = {\alpha}\nabla^2 u^{{*}}, .. _Eq:nonlin:splitting:RD:eq:decay: .. math:: \tag{599} \frac{\partial u^{{**}}}{\partial t} = f(u^{{**}}) {\thinspace .} Looking at these familiar problems, we may apply a :math:`\theta` rule (implicit) scheme for :ref:`(598) ` over one time step and avoid dealing with nonlinearities by applying an explicit scheme for :ref:`(599) ` over the same time step. Suppose we have some solution :math:`u` at time level :math:`t_n`. For flexibility, we define a :math:`\theta` method for the diffusion part :ref:`(598) ` by .. math:: [D_t u^{*} = {\alpha} (D_xD_x u^{*} + D_y D_y u^{*})]^{n+\theta}{\thinspace .} We use :math:`u^{n}` as initial condition for :math:`u^{*}`. 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 :math:`{\mathcal{O}(\Delta t^2)}` for Strang splitting, otherwise it is just :math:`{\mathcal{O}(\Delta t)}`. Higher-order methods for ODEs will therefore be a waste of work. The 2nd-order Adams-Bashforth method reads .. _Eq:_auto251: .. math:: \tag{600} u^{{**},n+1}_{i,j} = u^{{**},n}_{i,j} + \frac{1}{2}\Delta t\left( 3f(u^{{**}, n}_{i,j}, t_n) - f(u^{{**}, n-1}_{i,j}, t_{n-1}) \right) {\thinspace .} We can use a Forward Euler step to start the method, i.e, compute :math:`u^{{**},1}_{i,j}`. The algorithm goes like this: 1. Solve the diffusion problem for one time step as usual. 2. Solve the reaction ODEs at each mesh point in :math:`[t_n,t_n+\Delta t]`, using the diffusion solution in 1. as initial condition. The solution of the ODEs constitute 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 :math:`u'=f(u)`. This gives great flexibility in splitting methods. .. _nonlin:splitting:RD_linearR: 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. :math:`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 :math:`E = C h^r`, while differing in the step :math:`h` (being either :math:`dx^2` or :math:`dx`) and the convergence rate :math:`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: .. code-block:: python def convergence_rates(scheme='diffusion'): F = 0.5 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] for Nx in Nx_values: dx = L/Nx dt = F/a*dx**2 Nt = int(round(T/float(dt))) t = np.linspace(0, Nt*dt, Nt+1) # Mesh points, global time if scheme == 'diffusion': print 'Running 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 'Running 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 'Running 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) elif scheme == 'Strang_splitting_2ndOrder': print 'Running Strang splitting with 2nd order schemes...' 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) else: print 'Unknown scheme requested!' sys.exit(0) E.append(error) 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) The Fourier number is kept fixed throughout at :math:`F = 0.5`, being the stability limit with explicit schemes. Since :math:`\alpha` is considered a known constant and :math:`F = \frac{\alpha dt}{dx^2}`, :math:`dt` (or :math:`dx`) is readily computed for a given :math:`dx` (or :math:`dt`). The loop in ``convergence_rates`` runs over a chosen set of grid points implying a doubling of spatial resolution with each iteration. A solver ``diffusion_theta`` is used in each of the four solution approaches: .. code-block:: python 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``. It was implemented as .. code-block:: python 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 half, 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: .. code-block:: python 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 .. code-block:: python 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 .. code-block:: python 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 has second order convergence (:math:`r = 2`), while the remaining three approaches have first order convergence (:math:`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 :math:`f(u)=-\beta u` for a constant :math:`\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 :math:`x=0` and :math:`x=L`: .. math:: u = e^{-{\alpha} k^2 t - \beta t}\sin kx,\quad k = \frac{\pi}{L}{\thinspace .} This component fits the 1D PDE problem (:math:`f=0`). On complex form we can write .. math:: u = e^{-{\alpha} k^2 t - \beta t + ikx}, where :math:`i=\sqrt{-1}` and the imaginary part is taken as the physical solution. We refer to the section :ref:`diffu:pde1:analysis` and to the book [Ref02]_ for a discussion of exact numerical solutions to diffusion and decay problems, respectively. The key idea is to search for solutions :math:`A^ne^{ikx}` and determine :math:`A`. For the diffusion problem solved by a Forward Euler method one has .. math:: A = 1 - 4F\sin^p, where :math:`F={\alpha}\Delta t/\Delta x^2` is the mesh Fourier number and :math:`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 :math:`u'=-\beta u`, we have :math:`A=1 - q`, where :math:`q` is a dimensionless parameter reflecting the resolution in the decay problem: :math:`q = \beta\Delta t`. The original model problem can also be discretized by a Forward Euler scheme, .. math:: [D^+_t u = {\alpha} D_xD_x u - \beta u]^n_i{\thinspace .} Assuming :math:`A^ne^{ikx}` we find that .. math:: 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, .. math:: u^{n}_{i} = (1-4F\sin^2 p)u^{n-1}_i{\thinspace .} In the two stage algorithm, we first compute the diffusion step .. math:: u^{{*},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 .. math:: u^{{**},n+1} = (1-q)u^{{*},n+1} = (1-q)(1-4F\sin^2 p) u^{n-1}_i{\thinspace .} The splitting approximation over one step is therefore .. math:: E = 1 - 4F\sin^p -q - (1-q)(1-4F\sin^2 p) = -q(2 - F\sin^2 p))