{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Operator splitting methods\n", "
\n", "\n", "\n", "Operator splitting is a natural and old idea. When a PDE or system of PDEs\n", "contains different terms expressing different physics, it is natural to\n", "use different numerical methods for different physical processes. This can\n", "optimize and simplify the overall solution process. The idea was\n", "especially popularized in the context of the Navier-Stokes equations\n", "and reaction-diffusion PDEs. Common names for the technique are *operator\n", "splitting*, *fractional step* methods, and *split-step* methods. We shall\n", "stick to the former name.\n", "In the context of nonlinear differential equations, operator splitting\n", "can be used to isolate nonlinear terms and simplify the solution methods.\n", "\n", "A related technique, often known as dimensional splitting or\n", "alternating direction implicit (ADI) methods, is to split the spatial\n", "dimensions and solve a 2D or 3D problem as two or three\n", "consecutive 1D problems, but this type of splitting\n", "is not to be further considered here.\n", "\n", "## Ordinary operator splitting for ODEs\n", "
\n", "\n", "\n", "Consider first an ODE where the right-hand side is split into two\n", "terms:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "u' = f_0(u) + f_1(u)\\thinspace .\n", "\\label{_auto1} \\tag{1}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In case $f_0$ and $f_1$ are linear functions of $u$, $f_0=au$ and $f_1=bu$,\n", "we have $u(t)=Ie^{(a+b)t}$, if $u(0)=I$.\n", "When going one time step of length $\\Delta t$ from $t_n$ to $t_{n+1}$, we have" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u(t_{n+1}) = u(t_n)e^{(a+b)\\Delta t}\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This expression can be also be written as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u(t_{n+1}) = u(t_n)e^{a\\Delta t}e^{b\\Delta t},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "or" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "u^\\stepzero = u(t_n)e^{a\\Delta t},\n", "\\label{nonlin:splitting:ODE:step1} \\tag{2}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u(t_{n+1}) = u^\\stepzero e^{b\\Delta t}\n", "\\label{nonlin:splitting:ODE:step2} \\tag{3}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first step ([2](#nonlin:splitting:ODE:step1)) means solving\n", "$u'=f_0$ over a time interval $\\Delta t$ with $u(t_n)$ as start value.\n", "The second step ([3](#nonlin:splitting:ODE:step2)) means solving\n", "$u'=f_1$ over a time interval $\\Delta t$ with the value at the end\n", "of the first step as start value.\n", "That is, we progress the solution in two steps and solve\n", "two ODEs $u'=f_0$ and $u'=f_1$. The order of the equations is not\n", "important. From the derivation above we see that solving $u'=f_1$\n", "prior to $u'=f_0$ can equally well be done.\n", "\n", "The technique is exact if the ODEs are linear. For nonlinear ODEs it is\n", "only an approximate method with error $\\Delta t$. The technique can\n", "be extended to an arbitrary number of steps; i.e., we may split the PDE\n", "system into any number of subsystems. Examples will illuminate this principle.\n", "\n", "\n", "## Strang splitting for ODEs\n", "
\n", "\n", "\n", "The accuracy of the splitting method in the section [Ordinary operator splitting for ODEs](#nonlin:splitting:ODE)\n", "can be improved from $\\Oof{\\Delta t}$ to $\\Oof{\\Delta t^2}$ using so-called\n", "*Strang splitting*, where we take frac{1}{2} a step with the $f_0$ operator,\n", "a full step with the $f_1$ operator, and finally frac{1}{2} another step with\n", "the $f_0$ operator. During a time interval $\\Delta t$ the algorithm can\n", "be written as follows." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "\\frac{du^{\\stepzero}}{dt} &= f_0(u^{\\stepzero}),\n", "\\quad u^{\\stepzero}(t_n)=u(t_n),\n", "\\quad t\\in [t_n,t_n+\\frac{1}{2}\\Delta t],\\\\\n", "\\frac{du^{\\stephalf}}{dt} &= f_1(u^{\\stephalf}),\n", "\\quad u^{\\stephalf}(t_n)=u^{\\stepzero}(t_{n+\\frac{1}{2}}),\n", "\\quad t\\in [t_n,t_n+\\Delta t],\\\\\n", "\\frac{du^{\\stepone}}{dt} &= f_0(u^{\\stepone}),\n", "\\quad u^{\\stepone}(t_{n+\\frac{1}{2}})=u^{\\stephalf}(t_{n+1}),\n", "\\quad t\\in [t_n+\\frac{1}{2}\\Delta t, t_n+\\Delta t]\\thinspace .\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The global solution is set as $u(t_{n+1}) = u^{\\stepone}(t_{n+1})$.\n", "\n", "There is no use in combining higher-order methods with\n", "ordinary splitting since the error due to splitting is $\\Oof{\\Delta\n", "t}$, but for Strang splitting it makes sense to use schemes of order\n", "$\\Oof{\\Delta t^2}$.\n", "\n", "With the notation introduced for Strang splitting, we may express\n", "ordinary first-order splitting as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "\\frac{du^{\\stepzero}}{dt} &= f_0(u^{\\stepzero}),\\quad u^{\\stepzero}(t_n)=u(t_n),\n", "\\quad t\\in [t_n,t_n+\\Delta t],\\\\\n", "\\frac{du^{\\stepone}}{dt} &= f_1(u^{\\stepone}),\\quad u^{\\stepone}(t_n)=u^{\\stepzero}(t_{n+1}),\n", "\\quad t\\in [t_n,t_n+\\Delta t],\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "with global solution set as $u(t_{n+1}) = u^{\\stepone}(t_{n+1})$.\n", "\n", "## Example: Logistic growth\n", "
\n", "\n", "\n", "Let us split the (scaled) logistic equation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u'=u(1-u),\\quad u(0)=0.1,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "with solution $u=(9e^{-t}+1)^{-1}$, into" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u'=u - u^2 = f_0(u) + f_1(u), \\quad f_0(u)=u,\\quad f_1(u)=-u^2\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We solve $u'=f_0(u)$ and $u'=f_1(u)$ by a Forward Euler step.\n", "In addition, we add a method where we solve $u'=f_0(u)$ analytically,\n", "since the equation is actually $u'=u$ with solution $e^t$.\n", "The software that accompanies the following methods is the file\n", "[`split_logistic.py`](${src_nonlin}/split_logistic.py).\n", "\n", "### Splitting techniques\n", "\n", "Ordinary splitting takes a Forward Euler step for each of the ODEs\n", "according to" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{u^{\\stepzero,n+1} - u^{\\stepzero,n}}{\\Delta t} =\n", "f_0(u^{\\stepzero,n}),\\quad\n", "u^{\\stepzero,n} = u(t_n),\\quad t\\in [t_n,t_n+\\Delta t],\n", "\\label{_auto2} \\tag{4}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "\\frac{u^{\\stepone,n+1} - u^{\\stepone, n}}{\\Delta t} = f_1(u^{\\stepone,n}),\\quad\n", "u^{\\stepone,n} = u^{\\stepzero,n+1},\\quad t\\in [t_n,t_n+\\Delta t],\n", "\\label{_auto3} \\tag{5}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "with $u(t_{n+1}) = u^{\\stepone,n+1}$.\n", "\n", "Strang splitting takes the form" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{u^{\\stepzero,n+\\frac{1}{2}} - u^{\\stepzero,n}}{\\frac{1}{2}\\Delta t} =\n", "f_0(u^{\\stepzero,n}),\\quad\n", "u^{\\stepzero,n} = u(t_n),\\quad t\\in [t_n,t_n+\\frac{1}{2}\\Delta t],\n", "\\label{_auto4} \\tag{6}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "\\frac{u^{\\stephalf,n+1}-u^{\\stephalf,n}}{\\Delta t} =\n", "f_1(u^{\\stephalf,n}),\\quad\n", "u^{\\stephalf,n} = u^{\\stepzero, n+\\frac{1}{2}},\\quad t\\in [t_n,t_n+\\Delta t],\n", "\\label{_auto5} \\tag{7}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "\\frac{u^{\\stepone, n+1} - u^{\\stepone, n+\\frac{1}{2}}}{\\frac{1}{2}\\Delta t} =\n", "f_0(u^{\\stepone,n+\\frac{1}{2}}),\\quad\n", "u^{\\stepone,n+\\frac{1}{2}} = u^{\\stephalf,n+1},\\quad\n", "t\\in [t_n+\\frac{1}{2}\\Delta t, t_n+\\Delta t]\\thinspace .\n", "\\label{_auto6} \\tag{8}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Verbose implementation\n", "\n", "The following function computes four solutions arising from the Forward\n", "Euler method, ordinary splitting, Strang splitting, as well as Strang splitting\n", "with exact treatment of $u'=f_0(u)$:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "def solver(dt, T, f, f_0, f_1):\n", " \"\"\"\n", " Solve u'=f by the Forward Euler method and by ordinary and\n", " Strang splitting: f(u) = f_0(u) + f_1(u).\n", " \"\"\"\n", " Nt = int(round(T/float(dt)))\n", " t = np.linspace(0, Nt*dt, Nt+1)\n", " u_FE = np.zeros(len(t))\n", " u_split1 = np.zeros(len(t)) # 1st-order splitting\n", " u_split2 = np.zeros(len(t)) # 2nd-order splitting\n", " u_split3 = np.zeros(len(t)) # 2nd-order splitting w/exact f_0\n", "\n", " # Set initial values\n", " u_FE[0] = 0.1\n", " u_split1[0] = 0.1\n", " u_split2[0] = 0.1\n", " u_split3[0] = 0.1\n", "\n", " for n in range(len(t)-1):\n", " # Forward Euler method\n", " u_FE[n+1] = u_FE[n] + dt*f(u_FE[n])\n", "\n", " # --- Ordinary splitting ---\n", " # First step\n", " u_s_n = u_split1[n]\n", " u_s = u_s_n + dt*f_0(u_s_n)\n", " # Second step\n", " u_ss_n = u_s\n", " u_ss = u_ss_n + dt*f_1(u_ss_n)\n", " u_split1[n+1] = u_ss\n", "\n", " # --- Strang splitting ---\n", " # First step\n", " u_s_n = u_split2[n]\n", " u_s = u_s_n + dt/2.*f_0(u_s_n)\n", " # Second step\n", " u_sss_n = u_s\n", " u_sss = u_sss_n + dt*f_1(u_sss_n)\n", " # Third step\n", " u_ss_n = u_sss\n", " u_ss = u_ss_n + dt/2.*f_0(u_ss_n)\n", " u_split2[n+1] = u_ss\n", "\n", " # --- Strang splitting using exact integrator for u'=f_0 ---\n", " # First step\n", " u_s_n = u_split3[n]\n", " u_s = u_s_n*np.exp(dt/2.) # exact\n", " # Second step\n", " u_sss_n = u_s\n", " u_sss = u_sss_n + dt*f_1(u_sss_n)\n", " # Third step\n", " u_ss_n = u_sss\n", " u_ss = u_ss_n*np.exp(dt/2.) # exact\n", " u_split3[n+1] = u_ss\n", "\n", " return u_FE, u_split1, u_split2, u_split3, t" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Compact implementation\n", "\n", "We have used quite many lines for the steps in the splitting methods.\n", "Many will prefer to condense the code a bit, as done here:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Ordinary splitting\n", "u_s = u_split1[n] + dt*f_0(u_split1[n])\n", "u_split1[n+1] = u_s + dt*f_1(u_s)\n", "\n", "# Strang splitting\n", "u_s = u_split2[n] + dt/2.*f_0(u_split2[n])\n", "u_sss = u_s + dt*f_1(u_s)\n", "u_split2[n+1] = u_sss + dt/2.*f_0(u_sss)\n", "\n", "# Strang splitting using exact integrator for u'=f_0\n", "u_s = u_split3[n]*np.exp(dt/2.) # exact\n", "u_ss = u_s + dt*f_1(u_s)\n", "u_split3[n+1] = u_ss*np.exp(dt/2.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Results\n", "\n", "[Figure](#nonlin:splitting:ODE:logistic:fig) shows that the impact of\n", "splitting is significant. Interestingly, however,\n", "the Forward Euler method applied to the entire problem directly is much more\n", "accurate than any of the splitting schemes. We also see that Strang\n", "splitting is definitely more accurate than ordinary splitting and that\n", "it helps a bit to use an exact solution of $u'=f_0(u)$. With a large\n", "time step ($\\Delta t = 0.2$, left plot in [Figure](#nonlin:splitting:ODE:logistic:fig)), the asymptotic values are off\n", "by 20-30%. A more reasonable time step ($\\Delta t = 0.05$, right plot\n", "in [Figure](#nonlin:splitting:ODE:logistic:fig)) gives better\n", "results, but still the asymptotic values are up to 10% wrong.\n", "\n", "As technique for solving nonlinear ODEs, we realize that the present\n", "case study is not particularly promising, as the Forward Euler method\n", "both linearizes the original problem and provides a solution that is\n", "much more accurate than any of the splitting techniques.\n", "In complicated multi-physics settings, on the other hand, splitting\n", "may be the only feasible way to go, and sometimes you really need to apply\n", "different numerics to different parts of a PDE problem.\n", "But in very simple problems, like the logistic ODE, splitting is just an\n", "inferior technique. Still, the logistic ODE is ideal for introducing all\n", "the mathematical details and for investigating the behavior.\n", "\n", "\n", "\n", "
\n", "\n", "

Effect of ordinary and Strang splitting for the logistic equation.

\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Reaction-diffusion equation\n", "
\n", "\n", "Consider a diffusion equation coupled to chemical reactions modeled by\n", "a nonlinear term $f(u)$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial u}{\\partial t} = \\dfc\\nabla^2u + f(u)\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is a physical process composed of two individual processes:\n", "$u$ is the concentration of a substance that is locally generated by\n", "a chemical reaction $f(u)$, while $u$ is spreading in space because\n", "of diffusion. There are obviously two time scales: one for the chemical\n", "reaction and one for diffusion. Typically, fast chemical reactions require\n", "much finer time stepping than slower diffusion processes.\n", "mathcal{I}_t could therefore be advantageous to split\n", "the two physical effects in separate models and use different numerical methods\n", "for the two.\n", "\n", "A natural spitting in the present case is" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial u^{\\stepzero}}{\\partial t} = \\dfc\\nabla^2 u^{\\stepzero},\n", "\\label{nonlin:splitting:RD:eq:diffu} \\tag{9}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "\\frac{\\partial u^{\\stepone}}{\\partial t} = f(u^{\\stepone})\n", "\\label{nonlin:splitting:RD:eq:decay} \\tag{10}\\thinspace .\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Looking at these familiar problems,\n", "we may apply a $\\theta$ rule (implicit) scheme for\n", "([9](#nonlin:splitting:RD:eq:diffu)) over one time step\n", "and avoid dealing with nonlinearities\n", "by applying an explicit scheme for ([10](#nonlin:splitting:RD:eq:decay))\n", "over the same time step.\n", "\n", "Suppose we have some solution $u$ at time level $t_n$. For flexibility,\n", "we define a $\\theta$ method for the diffusion part\n", "([9](#nonlin:splitting:RD:eq:diffu)) by" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "[D_t u^\\stepzero =\n", "\\dfc (D_xD_x u^\\stepzero + D_y D_y u^\\stepzero)]^{n+\\theta}\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We use $u^{n}$ as initial condition for $u^\\stepzero$.\n", "\n", "The reaction part, which is defined at each mesh point (without coupling\n", "values in different mesh points), can employ any scheme for an ODE.\n", "Here we use an Adams-Bashforth method of order 2. Recall that the overall\n", "accuracy of the splitting method\n", "is maximum $\\Oof{\\Delta t^2}$ for Strang splitting,\n", "otherwise it is just $\\Oof{\\Delta t}$. Higher-order methods for ODEs will\n", "therefore be a waste of work. The 2nd-order Adams-Bashforth method reads" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "u^{\\stepone,n+1}_{i,j} = u^{\\stepone,n}_{i,j} +\n", "\\frac{1}{2}\\Delta t\\left( 3f(u^{\\stepone, n}_{i,j}, t_n) -\n", "f(u^{\\stepone, n-1}_{i,j}, t_{n-1})\n", "\\right)\n", "\\thinspace .\n", "\\label{_auto7} \\tag{11}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can use a Forward Euler step to start the method, i.e, compute\n", "$u^{\\stepone,1}_{i,j}$.\n", "\n", "The algorithm goes like this:\n", "\n", "1. Solve the diffusion problem for one time step as usual.\n", "\n", "2. Solve the reaction ODEs at each mesh point in $[t_n,t_n+\\Delta t]$,\n", " using the diffusion solution in 1. as initial condition.\n", " The solution of the ODEs constitutes the solution of the original problem\n", " at the end of each time step.\n", "\n", "We may use a much smaller time step when solving the reaction part, adapted\n", "to the dynamics of the problem $u'=f(u)$. This gives great flexibility in\n", "splitting methods.\n", "\n", "[hpl 1: Can you make an example? I have started in `src/split_diffu_react.py`.]\n", "\n", "## Example: Reaction-Diffusion with linear reaction term\n", "
\n", "\n", "\n", "The methods above may be explored in detail through a specific\n", "computational example in which we compute the convergence rates\n", "associated with four different solution approaches for the\n", "reaction-diffusion equation with a linear reaction term,\n", "i.e. $f(u)=-bu$. The methods comprise solving without splitting (just\n", "straight Forward Euler), ordinary splitting, first order Strang\n", "splitting, and second order Strang splitting. In all four methods, a\n", "standard centered difference approximation is used for the spatial\n", "second derivative. The methods share the error model $E = C h^r$,\n", "while differing in the step $h$ (being either $\\Delta x^2$ or $\\Delta\n", "x$) and the convergence rate $r$ (being either 1 or 2).\n", "\n", "All code commented below is found in the file\n", "[`split_diffu_react.py`](${src_nonlin}/split_diffu_react.py). When executed,\n", "a function `convergence_rates` is called, from which all convergence\n", "rate computations are handled:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def convergence_rates(scheme='diffusion'):\n", " \n", " F = 0.5 # Upper limit for FE (stability). For CN, this \n", " # limit does not apply, but for simplicity, we\n", " # choose F = 0.5 as the initial F value.\n", " T = 1.2\n", " a = 3.5\n", " b = 1\n", " L = 1.5\n", " k = np.pi/L\n", " \n", " def exact(x, t):\n", " '''exact sol. to: du/dt = a*d^2u/dx^2 - b*u'''\n", " return np.exp(-(a*k**2 + b)*t) * np.sin(k*x)\n", " \n", " def f(u, t):\n", " return -b*u \n", "\n", " def I(x):\n", " return exact(x, 0)\n", " \n", " global error # error computed in the user action function\n", " error = 0 \n", "\n", " # Convergence study\n", " def action(u, x, t, n):\n", " global error \n", " if n == 1: # New simulation, - reset error\n", " error = 0\n", " else:\n", " error = max(error, np.abs(u - exact(x, t[n])).max()) \n", "\n", " E = []\n", " h = []\n", " Nx_values = [10, 20, 40, 80] # i.e., dx halved each time\n", " for Nx in Nx_values: \n", " dx = L/Nx \n", " if scheme == 'Strang_splitting_2ndOrder':\n", " print 'Strang splitting with 2nd order schemes...'\n", " # In this case, E = C*h**r (with r = 2) and since \n", " # h = dx = K*dt, the ratio dt/dx must be constant.\n", " # To fulfill this demand, we must let F change\n", " # when dx changes. From F = a*dt/dx**2, it follows\n", " # that halving dx AND doubling F assures dt/dx const.\n", " # Initially, we simply choose F = 0.5.\n", "\n", " dt = F/a*dx**2\n", " #print 'dt/dx:', dt/dx \n", " Nt = int(round(T/float(dt)))\n", " t = np.linspace(0, Nt*dt, Nt+1) # global time \n", " Strang_splitting_2ndOrder(I=I, a=a, b=b, f=f, L=L, dt=dt, \n", " dt_Rfactor=1, F=F, t=t, T=T,\n", " user_action=action) \n", " h.append(dx)\n", " # prepare for next iteration (make F match dx/2)\n", " F = F*2 # assures dt/dx const. when dx = dx/2 \n", " else: \n", " # In these cases, E = C*h**r (with r = 1) and since \n", " # h = dx**2 = K*dt, the ratio dt/dx**2 must be constant.\n", " # This is fulfilled by choosing F = 0.5 (for FE stability)\n", " # and make sure that F, dx and dt comply to F = a*dt/dx**2. \n", " dt = F/a*dx**2\n", " Nt = int(round(T/float(dt)))\n", " t = np.linspace(0, Nt*dt, Nt+1) # global time\n", " if scheme == 'diffusion':\n", " print 'FE on whole eqn...'\n", " diffusion_theta(I, a, f, L, dt, F, t, T,\n", " step_no=0, theta=0,\n", " u_L=0, u_R=0, user_action=action) \n", " h.append(dx**2)\n", " elif scheme == 'ordinary_splitting':\n", " print 'Ordinary splitting...'\n", " ordinary_splitting(I=I, a=a, b=b, f=f, L=L, dt=dt, \n", " dt_Rfactor=1, F=F, t=t, T=T,\n", " user_action=action) \n", " h.append(dx**2)\n", " elif scheme == 'Strang_splitting_1stOrder':\n", " print 'Strang splitting with 1st order schemes...'\n", " Strang_splitting_1stOrder(I=I, a=a, b=b, f=f, L=L, dt=dt, \n", " dt_Rfactor=1, F=F, t=t, T=T,\n", " user_action=action) \n", " h.append(dx**2)\n", " else:\n", " print 'Unknown scheme requested!'\n", " sys.exit(0)\n", " \n", " #print 'dt/dx**2:', dt/dx**2 \n", " \n", " E.append(error)\n", " Nx *= 2 # Nx doubled gives dx/2\n", "\n", " print 'E:', E\n", " print 'h:', h\n", " \n", " # Convergence rates \n", " r = [np.log(E[i]/E[i-1])/np.log(h[i]/h[i-1])\n", " for i in range(1,len(Nx_values))]\n", " print 'Computed rates:', r\n", "\n", "\n", "if __name__ == '__main__':\n", " \n", " schemes = ['diffusion', \n", " 'ordinary_splitting', \n", " 'Strang_splitting_1stOrder',\n", " 'Strang_splitting_2ndOrder']\n", " \n", " for scheme in schemes:\n", " convergence_rates(scheme=scheme)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, with respect to the error ($E = C h^r$), the Forward Euler\n", "scheme, the ordinary splitting scheme and first order Strang splitting\n", "scheme are all first order ($r = 1$), with a step $h = \\Delta x^2 =\n", "K^{-1}\\Delta t$, where $K$ is some constant. This implies that the\n", "*ratio* $\\frac{\\Delta t}{\\Delta x^2}$ must be held constant during\n", "convergence rate calculations. Furthermore, the Fourier number $F =\n", "\\frac{\\alpha \\Delta t}{\\Delta x^2}$ is upwards limited to $F = 0.5$,\n", "being the stability limit with explicit schemes. Thus, in these cases,\n", "we use the fixed value of $F$ and a given (but changing) spatial\n", "resolution $\\Delta x$ to compute the corresponding value of $\\Delta t$\n", "according to the expression for $F$. This assures that $\\frac{\\Delta\n", "t}{\\Delta x^2}$ is kept constant. The loop in `convergence_rates` runs\n", "over a chosen set of grid points (`Nx_values`) which gives a doubling\n", "of spatial resolution with each iteration ($\\Delta x$ is halved).\n", "\n", "For the second order Strang splitting scheme, we have $r = 2$ and a\n", "step $h = \\Delta x = K^{-1}\\Delta t$, where $K$ again is some\n", "constant. In this case, it is thus the ratio $\\frac{\\Delta t}{\\Delta\n", "x}$ that must be held constant during the convergence rate\n", "calculations. From the expression for $F$, it is clear then that $F$\n", "must change with each halving of $\\Delta x$. In fact, if $F$ is\n", "doubled each time $\\Delta x$ is halved, the ratio $\\frac{\\Delta\n", "t}{\\Delta x}$ will be constant (this follows, e.g., from the\n", "expression for $F$). This is utilized in our code.\n", "\n", "A solver `diffusion_theta` is used in each of the four solution approaches:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def diffusion_theta(I, a, f, L, dt, F, t, T, step_no, theta=0.5, \n", " u_L=0, u_R=0, user_action=None):\n", " \"\"\"\n", " Full solver for the model problem using the theta-rule\n", " difference approximation in time (no restriction on F,\n", " i.e., the time step when theta >= 0.5). Vectorized \n", " implementation and sparse (tridiagonal) coefficient matrix.\n", " Note that t always covers the whole global time interval, whether\n", " splitting is the case or not. T, on the other hand, is \n", " the end of the global time interval if there is no split,\n", " but if splitting, we use T=dt. When splitting, step_no \n", " keeps track of the time step number (for lookup in t). \n", " \"\"\"\n", "\n", " Nt = int(round(T/float(dt)))\n", " dx = np.sqrt(a*dt/F) \n", " Nx = int(round(L/dx))\n", " x = np.linspace(0, L, Nx+1) # Mesh points in space\n", " # Make sure dx and dt are compatible with x and t\n", " dx = x[1] - x[0]\n", " dt = t[1] - t[0]\n", "\n", " u = np.zeros(Nx+1) # solution array at t[n+1]\n", " u_1 = np.zeros(Nx+1) # solution at t[n]\n", "\n", " # Representation of sparse matrix and right-hand side\n", " diagonal = np.zeros(Nx+1)\n", " lower = np.zeros(Nx)\n", " upper = np.zeros(Nx)\n", " b = np.zeros(Nx+1)\n", "\n", " # Precompute sparse matrix (scipy format)\n", " Fl = F*theta\n", " Fr = F*(1-theta)\n", " diagonal[:] = 1 + 2*Fl\n", " lower[:] = -Fl #1\n", " upper[:] = -Fl #1\n", " # Insert boundary conditions\n", " diagonal[0] = 1\n", " upper[0] = 0\n", " diagonal[Nx] = 1\n", " lower[-1] = 0\n", "\n", " diags = [0, -1, 1]\n", " A = scipy.sparse.diags(\n", " diagonals=[diagonal, lower, upper],\n", " offsets=[0, -1, 1], shape=(Nx+1, Nx+1),\n", " format='csr')\n", " #print A.todense()\n", "\n", " # Allow f to be None or 0\n", " if f is None or f == 0:\n", " f = lambda x, t: np.zeros((x.size)) \\\n", " if isinstance(x, np.ndarray) else 0\n", "\n", " # Set initial condition \n", " if isinstance(I, np.ndarray): # I is an array\n", " u_1 = np.copy(I)\n", " else: # I is a function\n", " for i in range(0, Nx+1):\n", " u_1[i] = I(x[i])\n", "\n", " if user_action is not None:\n", " user_action(u_1, x, t, step_no+0) \n", "\n", " # Time loop\n", " for n in range(0, Nt):\n", " b[1:-1] = u_1[1:-1] + \\\n", " Fr*(u_1[:-2] - 2*u_1[1:-1] + u_1[2:]) + \\\n", " dt*theta*f(u_1[1:-1], t[step_no+n+1]) + \\\n", " dt*(1-theta)*f(u_1[1:-1], t[step_no+n])\n", " b[0] = u_L; b[-1] = u_R # boundary conditions\n", " u[:] = scipy.sparse.linalg.spsolve(A, b)\n", "\n", " if user_action is not None:\n", " user_action(u, x, t, step_no+(n+1))\n", "\n", " # Update u_1 before next step\n", " u_1, u = u, u_1\n", "\n", " # u is now contained in u_1 (swapping)\n", " return u_1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the no splitting approach with Forward Euler in time, this solver\n", "handles both the diffusion and the reaction term. When splitting,\n", "`diffusion_theta` takes care of the diffusion term only, while the\n", "reaction term is handled either by a Forward Euler scheme in\n", "`reaction_FE`, or by a second order Adams-Bashforth scheme from\n", "Odespy. The `reaction_FE` function covers one complete time step `dt`\n", "during ordinary splitting, while Strang splitting (both first and\n", "second order) applies it with `dt/2` twice during each time step\n", "`dt`. Since the reaction term typically represents a much faster\n", "process than the diffusion term, a further refinement of the time step\n", "is made possible in `reaction_FE`. mathcal{I}_t was implemented as" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def reaction_FE(I, f, L, Nx, dt, dt_Rfactor, t, step_no, \n", " user_action=None):\n", " \"\"\"Reaction solver, Forward Euler method.\n", " Note the at t covers the whole global time interval. \n", " dt is either one complete,or one frac{1}{2}, of the step in the \n", " diffusion part, i.e. there is a local time interval \n", " [0, dt] or [0, dt/2] that the reaction_FE\n", " deals with each time it is called. step_no keeps\n", " track of the (global) time step number (required \n", " for lookup in t). \n", " \"\"\"\n", "\n", " u = np.copy(I) \n", " dt_local = dt/float(dt_Rfactor) \n", " Nt_local = int(round(dt/float(dt_local))) \n", " x = np.linspace(0, L, Nx+1) \n", " \n", " for n in range(Nt_local):\n", " time = t[step_no] + n*dt_local\n", " u[1:Nx] = u[1:Nx] + dt_local*f(u[1:Nx], time) \n", " \n", " # BC already inserted in diffusion step, i.e. no action here \n", " return u" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With the ordinary splitting approach, each time step `dt` is covered\n", "twice. First computing the impact of the reaction term, then the\n", "contribution from the diffusion term:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def ordinary_splitting(I, a, b, f, L, dt, \n", " dt_Rfactor, F, t, T, \n", " user_action=None):\n", " '''1st order scheme, i.e. Forward Euler is enough for both\n", " the diffusion and the reaction part. The time step dt is \n", " given for the diffusion step, while the time step for the \n", " reaction part is found as dt/dt_Rfactor, where dt_Rfactor >= 1.\n", " '''\n", " Nt = int(round(T/float(dt)))\n", " dx = np.sqrt(a*dt/F)\n", " Nx = int(round(L/dx))\n", " x = np.linspace(0, L, Nx+1) # Mesh points in space\n", " u = np.zeros(Nx+1)\n", " \n", " # Set initial condition u(x,0) = I(x)\n", " for i in range(0, Nx+1):\n", " u[i] = I(x[i])\n", "\n", " # In the following loop, each time step is \"covered twice\", \n", " # first for reaction, then for diffusion\n", " for n in range(0, Nt): \n", " # Reaction step (potentially many smaller steps within dt)\n", " u_s = reaction_FE(I=u, f=f, L=L, Nx=Nx, \n", " dt=dt, dt_Rfactor=dt_Rfactor, \n", " t=t, step_no=n, \n", " user_action=None) \n", " \n", " u = diffusion_theta(I=u_s, a=a, f=0, L=L, dt=dt, F=F,\n", " t=t, T=dt, step_no=n, theta=0,\n", " u_L=0, u_R=0, user_action=None) \n", " \n", " if user_action is not None:\n", " user_action(u, x, t, n+1)\n", "\n", " return" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the two Strang splitting approaches, each time step `dt` is\n", "handled by first computing the reaction step for (the first) `dt/2`,\n", "followed by a diffusion step `dt`, before the reaction step is treated\n", "once again for (the remaining) `dt/2`. Since first order Strang\n", "splitting is no better than first order accurate, both the reaction\n", "and diffusion steps are computed explicitly. The solver was\n", "implemented as" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "def Strang_splitting_1stOrder(I, a, b, f, L, dt, dt_Rfactor, \n", " F, t, T, user_action=None):\n", " '''Strang splitting while still using FE for the reaction\n", " step and for the diffusion step. Gives 1st order scheme.\n", " The time step dt is given for the diffusion step, while \n", " the time step for the reaction part is found as \n", " 0.5*dt/dt_Rfactor, where dt_Rfactor >= 1. Introduce an \n", " extra time mesh t2 for the reaction part, since it steps dt/2.\n", " '''\n", " Nt = int(round(T/float(dt)))\n", " t2 = np.linspace(0, Nt*dt, (Nt+1)+Nt) # Mesh points in diff \n", " dx = np.sqrt(a*dt/F)\n", " Nx = int(round(L/dx))\n", " x = np.linspace(0, L, Nx+1) \n", " u = np.zeros(Nx+1)\n", " \n", " # Set initial condition u(x,0) = I(x)\n", " for i in range(0, Nx+1):\n", " u[i] = I(x[i])\n", "\n", " for n in range(0, Nt): \n", " # Reaction step (1/2 dt: from t_n to t_n+1/2)\n", " # (potentially many smaller steps within dt/2)\n", " u_s = reaction_FE(I=u, f=f, L=L, Nx=Nx, \n", " dt=dt/2.0, dt_Rfactor=dt_Rfactor, \n", " t=t2, step_no=2*n, \n", " user_action=None)\n", " # Diffusion step (1 dt: from t_n to t_n+1)\n", " u_sss = diffusion_theta(I=u_s, a=a, f=0, L=L, dt=dt, F=F,\n", " t=t, T=dt, step_no=n, theta=0,\n", " u_L=0, u_R=0, user_action=None) \n", " # Reaction step (1/2 dt: from t_n+1/2 to t_n+1)\n", " # (potentially many smaller steps within dt/2)\n", " u = reaction_FE(I=u_sss, f=f, L=L, Nx=Nx, \n", " dt=dt/2.0, dt_Rfactor=dt_Rfactor, \n", " t=t2, step_no=2*n+1, \n", " user_action=None)\n", " \n", " if user_action is not None:\n", " user_action(u, x, t, n+1)\n", "\n", " return" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The second order version of the Strang splitting approach utilizes a\n", "second order Adams-Bashforth solver for the reaction part and a\n", "Crank-Nicolson scheme for the diffusion part. The solver has the same\n", "structure as the one for first order Strang splitting and was\n", "implemented as" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "def Strang_splitting_2ndOrder(I, a, b, f, L, dt, dt_Rfactor, \n", " F, t, T, user_action=None):\n", " '''Strang splitting using Crank-Nicolson for the diffusion\n", " step (theta-rule) and Adams-Bashforth 2 for the reaction step. \n", " Gives 2nd order scheme. Introduce an extra time mesh t2 for \n", " the reaction part, since it steps dt/2.\n", " '''\n", " import odespy\n", " Nt = int(round(T/float(dt)))\n", " t2 = np.linspace(0, Nt*dt, (Nt+1)+Nt) # Mesh points in diff \n", " dx = np.sqrt(a*dt/F)\n", " Nx = int(round(L/dx))\n", " x = np.linspace(0, L, Nx+1) \n", " u = np.zeros(Nx+1)\n", " \n", " # Set initial condition u(x,0) = I(x)\n", " for i in range(0, Nx+1):\n", " u[i] = I(x[i])\n", "\n", " reaction_solver = odespy.AdamsBashforth2(f) \n", "\n", " for n in range(0, Nt): \n", " # Reaction step (1/2 dt: from t_n to t_n+1/2)\n", " # (potentially many smaller steps within dt/2)\n", " reaction_solver.set_initial_condition(u)\n", " t_points = np.linspace(0, dt/2.0, dt_Rfactor+1)\n", " u_AB2, t_ = reaction_solver.solve(t_points) # t_ not needed\n", " u_s = u_AB2[-1,:] # pick sol at last point in time\n", "\n", " # Diffusion step (1 dt: from t_n to t_n+1)\n", " u_sss = diffusion_theta(I=u_s, a=a, f=0, L=L, dt=dt, F=F,\n", " t=t, T=dt, step_no=n, theta=0.5,\n", " u_L=0, u_R=0, user_action=None) \n", " # Reaction step (1/2 dt: from t_n+1/2 to t_n+1)\n", " # (potentially many smaller steps within dt/2)\n", " reaction_solver.set_initial_condition(u_sss)\n", " t_points = np.linspace(0, dt/2.0, dt_Rfactor+1)\n", " u_AB2, t_ = reaction_solver.solve(t_points) # t_ not needed\n", " u = u_AB2[-1,:] # pick sol at last point in time\n", " \n", " if user_action is not None:\n", " user_action(u, x, t, n+1)\n", "\n", " return" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When executing `split_diffu_react.py`, we find that the estimated\n", "convergence rates are as expected. The second order Strang splitting\n", "gives the least error (about $4e^{-5}$) and has second order\n", "convergence ($r = 2$), while the remaining three approaches have first\n", "order convergence ($r = 1$).\n", "\n", "\n", "\n", "## Analysis of the splitting method\n", "\n", "Let us address a linear PDE problem for which we can develop analytical\n", "solutions of the discrete equations, with and without splitting, and discuss\n", "these. Choosing $f(u)=-\\beta u$ for a constant $\\beta$ gives a linear\n", "problem. We use the Forward Euler method for both the PDE and ODE problems.\n", "\n", "We seek a 1D Fourier wave component solution of the problem, assuming\n", "homogeneous Dirichlet conditions at $x=0$ and $x=L$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u = e^{-\\dfc k^2 t - \\beta t}\\sin kx,\\quad k = \\frac{\\pi}{L}\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This component fits the 1D PDE problem ($f=0$). On complex form we can\n", "write" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u = e^{-\\dfc k^2 t - \\beta t + ikx},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $i=\\sqrt{-1}$ and the imaginary part is taken as the physical solution.\n", "\n", "We refer to the section \"Analysis of schemes for\n", "the diffusion equation\": \"\" in [[Langtangen_Linge]](#Langtangen_Linge) and to\n", "the book [[Langtangen_decay]](#Langtangen_decay) for a discussion of exact numerical\n", "solutions to diffusion and decay problems, respectively. The key idea\n", "is to search for solutions $A^ne^{ikx}$ and determine $A$. For the\n", "diffusion problem solved by a Forward Euler method one has" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "A = 1 - 4F\\sin^p,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $F=\\dfc\\Delta t/\\Delta x^2$ is the mesh Fourier number and $p=k\\Delta x/2$\n", "is a dimensionless number reflecting the spatial resolution (number of points\n", "per wave length in space). For the decay problem $u'=-\\beta u$, we have\n", "$A=1 - q$, where $q$ is a dimensionless parameter for the resolution\n", "in the decay problem: $q = \\beta\\Delta t$.\n", "\n", "The original model problem can also be discretized by a Forward Euler scheme," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "[D^+_t u = \\dfc D_xD_x u - \\beta u]^n_i\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Assuming $A^ne^{ikx}$ we find that" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u^n_i = (1 - 4F\\sin^p -q)^n\\sin kx\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We are particularly interested in what happens at one time step. That is," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u^{n}_{i} = (1-4F\\sin^2 p)u^{n-1}_i\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the two stage algorithm, we first compute the diffusion step" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u^{\\stepzero,n+1}_i = (1 - 4F\\sin^2 p)u^{n-1}_i\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we use this as input to the decay algorithm and arrive at" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u^{\\stepone,n+1} = (1-q)u^{\\stepzero,n+1} = (1-q)(1-4F\\sin^2 p) u^{n-1}_i\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The splitting approximation over one step is therefore" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E = 1 - 4F\\sin^p -q - (1-q)(1-4F\\sin^2 p) = -q(2 - F\\sin^2 p))\\thinspace .\n", "$$" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.3" } }, "nbformat": 4, "nbformat_minor": 4 }