.. !split .. _trunc:decay: Exponential decay ODEs ====================== .. index:: decay ODE We shall now compute the truncation error of a finite difference scheme for a differential equation. Our first problem involves the following the linear ODE modeling exponential decay, .. _Eq:trunc:decay:ode: .. math:: \tag{680} u'(t)=-au(t){\thinspace .} .. _trunc:decay:FE: Forward Euler scheme (2) --------------------------------- We begin with the Forward Euler scheme for discretizing :ref:`(680) `: .. _Eq:trunc:decay:FE:scheme: .. math:: \tag{681} \lbrack D_t^+ u = -au \rbrack^n {\thinspace .} The idea behind the truncation error computation is to insert the exact solution :math:`{u_{\small\mbox{e}}}` of the differential equation problem :ref:`(680) ` in the discrete equations :ref:`(681) ` and find the residual that arises because :math:`{u_{\small\mbox{e}}}` does not solve the discrete equations. Instead, :math:`{u_{\small\mbox{e}}}` solves the discrete equations with a residual :math:`R^n`: .. _Eq:trunc:decay:FE:uex: .. math:: \tag{682} [D_t^+ {u_{\small\mbox{e}}} + a{u_{\small\mbox{e}}} = R]^n {\thinspace .} From :ref:`(664) `-:ref:`(665) ` it follows that .. math:: [D_t^+ {u_{\small\mbox{e}}}]^n = {u_{\small\mbox{e}}}'(t_n) + \frac{1}{2}{u_{\small\mbox{e}}}''(t_n)\Delta t + {\mathcal{O}(\Delta t^2)}, which inserted in :ref:`(682) ` results in .. math:: {u_{\small\mbox{e}}}'(t_n) + \frac{1}{2}{u_{\small\mbox{e}}}''(t_n)\Delta t + {\mathcal{O}(\Delta t^2)} + a{u_{\small\mbox{e}}}(t_n) = R^n {\thinspace .} Now, :math:`{u_{\small\mbox{e}}}'(t_n) + a{u_{\small\mbox{e}}}^n = 0` since :math:`{u_{\small\mbox{e}}}` solves the differential equation. The remaining terms constitute the residual: .. _Eq:trunc:decay:FE:R: .. math:: \tag{683} R^n = \frac{1}{2}{u_{\small\mbox{e}}}''(t_n)\Delta t + {\mathcal{O}(\Delta t^2)} {\thinspace .} This is the truncation error :math:`R^n` of the Forward Euler scheme. Because :math:`R^n` is proportional to :math:`\Delta t`, we say that the Forward Euler scheme is of first order in :math:`\Delta t`. However, the truncation error is just one error measure, and it is not equal to the true error :math:`{u_{\small\mbox{e}}}^n - u^n`. For this simple model problem we can compute a range of different error measures for the Forward Euler scheme, including the true error :math:`{u_{\small\mbox{e}}}^n - u^n`, and all of them have dominating terms proportional to :math:`\Delta t`. .. _trunc:decay:CN: Crank-Nicolson scheme (2) ---------------------------------- For the Crank-Nicolson scheme, .. _Eq:trunc:decay:CN:scheme: .. math:: \tag{684} [D_t u = -au]^{n+\frac{1}{2}}, we compute the truncation error by inserting the exact solution of the ODE and adding a residual :math:`R`, .. _Eq:trunc:decay:CN:scheme:R: .. math:: \tag{685} [D_t {u_{\small\mbox{e}}} + a\overline{{u_{\small\mbox{e}}}}^{t} = R]^{n+\frac{1}{2}} {\thinspace .} The term :math:`[D_t{u_{\small\mbox{e}}}]^{n+\frac{1}{2}}` is easily computed from :ref:`(658) `-:ref:`(659) ` by replacing :math:`n` with :math:`n+{\frac{1}{2}}` in the formula, .. math:: \lbrack D_t{u_{\small\mbox{e}}}\rbrack^{n+\frac{1}{2}} = {u_{\small\mbox{e}}}'(t_{n+\frac{1}{2}}) + \frac{1}{24}{u_{\small\mbox{e}}}'''(t_{n+\frac{1}{2}})\Delta t^2 + {\mathcal{O}(\Delta t^4)}{\thinspace .} The arithmetic mean is related to :math:`u(t_{n+\frac{1}{2}})` by :ref:`(674) `-:ref:`(675) ` so .. math:: [a\overline{{u_{\small\mbox{e}}}}^{t}]^{n+\frac{1}{2}} = {u_{\small\mbox{e}}}(t_{n+\frac{1}{2}}) + \frac{1}{8}{u_{\small\mbox{e}}}''(t_{n})\Delta t^2 + + {\mathcal{O}(\Delta t^4)}{\thinspace .} Inserting these expressions in :ref:`(685) ` and observing that :math:`{u_{\small\mbox{e}}}'(t_{n+\frac{1}{2}}) +a{u_{\small\mbox{e}}}^{n+\frac{1}{2}} = 0`, because :math:`{u_{\small\mbox{e}}}(t)` solves the ODE :math:`u'(t)=-au(t)` at any point :math:`t`, we find that .. _Eq:_auto256: .. math:: \tag{686} R^{n+\frac{1}{2}} = \left( \frac{1}{24}{u_{\small\mbox{e}}}'''(t_{n+\frac{1}{2}}) + \frac{1}{8}{u_{\small\mbox{e}}}''(t_{n}) \right)\Delta t^2 + {\mathcal{O}(\Delta t^4)} Here, the truncation error is of second order because the leading term in :math:`R` is proportional to :math:`\Delta t^2`. At this point it is wise to redo some of the computations above to establish the truncation error of the Backward Euler scheme, see :ref:`trunc:exer:decay:BE`. .. _trunc:decay:theta: The :math:`\theta`-rule ----------------------- We may also compute the truncation error of the :math:`\theta`-rule, .. math:: [\bar D_t u = -a\overline{u}^{t,\theta}]^{n+\theta} {\thinspace .} Our computational task is to find :math:`R^{n+\theta}` in .. math:: [\bar D_t {u_{\small\mbox{e}}} + a\overline{{u_{\small\mbox{e}}}}^{t,\theta} = R]^{n+\theta} {\thinspace .} From :ref:`(666) `-:ref:`(667) ` and :ref:`(672) `-:ref:`(673) ` we get expressions for the terms with :math:`{u_{\small\mbox{e}}}`. Using that :math:`{u_{\small\mbox{e}}}'(t_{n+\theta}) + a{u_{\small\mbox{e}}}(t_{n+\theta})=0`, we end up with .. math:: R^{n+\theta} = ({\frac{1}{2}}-\theta){u_{\small\mbox{e}}}''(t_{n+\theta})\Delta t + \frac{1}{2}\theta (1-\theta){u_{\small\mbox{e}}}''(t_{n+\theta})\Delta t^2 + \nonumber .. _Eq:_auto257: .. math:: \tag{687} \frac{1}{2}(\theta^2 -\theta + 3){u_{\small\mbox{e}}}'''(t_{n+\theta})\Delta t^2 + {\mathcal{O}(\Delta t^3)} For :math:`\theta =\frac{1}{2}` the first-order term vanishes and the scheme is of second order, while for :math:`\theta\neq \frac{1}{2}` we only have a first-order scheme. .. _trunc:decay:software: Using symbolic software ----------------------- The previously mentioned ``truncation_error`` module can be used to automate the Taylor series expansions and the process of collecting terms. Here is an example on possible use: .. code-block:: python from truncation_error import DiffOp from sympy import * def decay(): u, a = symbols('u a') diffop = DiffOp(u, independent_variable='t', num_terms_Taylor_series=3) D1u = diffop.D(1) # symbol for du/dt ODE = D1u + a*u # define ODE # Define schemes FE = diffop['Dtp'] + a*u CN = diffop['Dt' ] + a*u BE = diffop['Dtm'] + a*u theta = diffop['barDt'] + a*diffop['weighted_arithmetic_mean'] theta = sm.simplify(sm.expand(theta)) # Residuals (truncation errors) R = {'FE': FE-ODE, 'BE': BE-ODE, 'CN': CN-ODE, 'theta': theta-ODE} return R The returned dictionary becomes .. code-block:: text decay: { 'BE': D2u*dt/2 + D3u*dt**2/6, 'FE': -D2u*dt/2 + D3u*dt**2/6, 'CN': D3u*dt**2/24, 'theta': -D2u*a*dt**2*theta**2/2 + D2u*a*dt**2*theta/2 - D2u*dt*theta + D2u*dt/2 + D3u*a*dt**3*theta**3/3 - D3u*a*dt**3*theta**2/2 + D3u*a*dt**3*theta/6 + D3u*dt**2*theta**2/2 - D3u*dt**2*theta/2 + D3u*dt**2/6, } The results are in correspondence with our hand-derived expressions. .. _trunc:decay:estimate:R: Empirical verification of the truncation error ---------------------------------------------- The task of this section is to demonstrate how we can compute the truncation error :math:`R` numerically. For example, the truncation error of the Forward Euler scheme applied to the decay ODE :math:`u'=-ua` is .. _Eq:trunc:decay:FE:R:comp: .. math:: \tag{688} R^n = [D_t^+{u_{\small\mbox{e}}} + a{u_{\small\mbox{e}}}]^n {\thinspace .} If we happen to know the exact solution :math:`{u_{\small\mbox{e}}}(t)`, we can easily evaluate :math:`R^n` from the above formula. To estimate how :math:`R` varies with the discretization parameter :math:`\Delta t`, which has been our focus in the previous mathematical derivations, we first make the assumption that :math:`R=C\Delta t^r` for appropriate constants :math:`C` and :math:`r` and small enough :math:`\Delta t`. The rate :math:`r` can be estimated from a series of experiments where :math:`\Delta t` is varied. Suppose we have :math:`m` experiments :math:`(\Delta t_i, R_i)`, :math:`i=0,\ldots,m-1`. For two consecutive experiments :math:`(\Delta t_{i-1}, R_{i-1})` and :math:`(\Delta t_i, R_i)`, a corresponding :math:`r_{i-1}` can be estimated by .. _Eq:trunc:R:empir1: .. math:: \tag{689} r_{i-1} = \frac{\ln (R_{i-1}/R_i)}{\ln (\Delta t_{i-1}/\Delta t_i)}, for :math:`i=1,\ldots,m-1`. Note that the truncation error :math:`R_i` varies through the mesh, so :ref:`(689) ` is to be applied pointwise. A complicating issue is that :math:`R_i` and :math:`R_{i-1}` refer to different meshes. Pointwise comparisons of the truncation error at a certain point in all meshes therefore requires any computed :math:`R` to be restricted to the *coarsest mesh* and that all finer meshes contain all the points in the coarsest mesh. Suppose we have :math:`N_0` intervals in the coarsest mesh. Inserting a superscript :math:`n` in :ref:`(689) `, where :math:`n` counts mesh points in the coarsest mesh, :math:`n=0,\ldots,N_0`, leads to the formula .. _Eq:trunc:R:empir2: .. math:: \tag{690} r_{i-1}^n = \frac{\ln (R_{i-1}^n/R_i^n)}{\ln (\Delta t_{i-1}/\Delta t_i)} {\thinspace .} Experiments are most conveniently defined by :math:`N_0` and a number of refinements :math:`m`. Suppose each mesh has twice as many cells :math:`N_i` as the previous one: .. math:: N_i = 2^iN_0,\quad \Delta t_i = TN_i^{-1}, where :math:`[0,T]` is the total time interval for the computations. Suppose the computed :math:`R_i` values on the mesh with :math:`N_i` intervals are stored in an array ``R[i]`` (``R`` being a list of arrays, one for each mesh). Restricting this :math:`R_i` function to the coarsest mesh means extracting every :math:`N_i/N_0` point and is done as follows: .. code-block:: python stride = N[i]/N_0 R[i] = R[i][::stride] The quantity ``R[i][n]`` now corresponds to :math:`R_i^n`. In addition to estimating :math:`r` for the pointwise values of :math:`R=C\Delta t^r`, we may also consider an integrated quantity on mesh :math:`i`, .. _Eq:_auto258: .. math:: \tag{691} R_{I,i} = \left(\Delta t_i\sum_{n=0}^{N_i} (R_i^n)^2\right)^\frac{1}{2}\approx \int_0^T R_i(t)dt {\thinspace .} The sequence :math:`R_{I,i}`, :math:`i=0,\ldots,m-1`, is also expected to behave as :math:`C\Delta t^r`, with the same :math:`r` as for the pointwise quantity :math:`R`, as :math:`\Delta t\rightarrow 0`. The function below computes the :math:`R_i` and :math:`R_{I,i}` quantities, plots them and compares with the theoretically derived truncation error (``R_a``) if available. .. code-block:: python import numpy as np import scitools.std as plt def estimate(truncation_error, T, N_0, m, makeplot=True): """ Compute the truncation error in a problem with one independent variable, using m meshes, and estimate the convergence rate of the truncation error. The user-supplied function truncation_error(dt, N) computes the truncation error on a uniform mesh with N intervals of length dt:: R, t, R_a = truncation_error(dt, N) where R holds the truncation error at points in the array t, and R_a are the corresponding theoretical truncation error values (None if not available). The truncation_error function is run on a series of meshes with 2**i*N_0 intervals, i=0,1,...,m-1. The values of R and R_a are restricted to the coarsest mesh. and based on these data, the convergence rate of R (pointwise) and time-integrated R can be estimated empirically. """ N = [2**i*N_0 for i in range(m)] R_I = np.zeros(m) # time-integrated R values on various meshes R = [None]*m # time series of R restricted to coarsest mesh R_a = [None]*m # time series of R_a restricted to coarsest mesh dt = np.zeros(m) legends_R = []; legends_R_a = [] # all legends of curves for i in range(m): dt[i] = T/float(N[i]) R[i], t, R_a[i] = truncation_error(dt[i], N[i]) R_I[i] = np.sqrt(dt[i]*np.sum(R[i]**2)) if i == 0: t_coarse = t # the coarsest mesh stride = N[i]/N_0 R[i] = R[i][::stride] # restrict to coarsest mesh R_a[i] = R_a[i][::stride] if makeplot: plt.figure(1) plt.plot(t_coarse, R[i], log='y') legends_R.append('N=%d' % N[i]) plt.hold('on') plt.figure(2) plt.plot(t_coarse, R_a[i] - R[i], log='y') plt.hold('on') legends_R_a.append('N=%d' % N[i]) if makeplot: plt.figure(1) plt.xlabel('time') plt.ylabel('pointwise truncation error') plt.legend(legends_R) plt.savefig('R_series.png') plt.savefig('R_series.pdf') plt.figure(2) plt.xlabel('time') plt.ylabel('pointwise error in estimated truncation error') plt.legend(legends_R_a) plt.savefig('R_error.png') plt.savefig('R_error.pdf') # Convergence rates r_R_I = convergence_rates(dt, R_I) print 'R integrated in time; r:', print ' '.join(['%.1f' % r for r in r_R_I]) R = np.array(R) # two-dim. numpy array r_R = [convergence_rates(dt, R[:,n])[-1] for n in range(len(t_coarse))] The first ``makeplot`` block demonstrates how to build up two figures in parallel, using ``plt.figure(i)`` to create and switch to figure number ``i.`` Figure numbers start at 1. A logarithmic scale is used on the :math:`y` axis since we expect that :math:`R` as a function of time (or mesh points) is exponential. The reason is that the theoretical estimate :ref:`(683) ` contains :math:`{u_{\small\mbox{e}}}''`, which for the present model goes like :math:`e^{-at}`. Taking the logarithm makes a straight line. The code follows closely the previously stated mathematical formulas, but the statements for computing the convergence rates might deserve an explanation. The generic help function ``convergence_rate(h, E)`` computes and returns :math:`r_{i-1}`, :math:`i=1,\ldots,m-1` from :ref:`(690) `, given :math:`\Delta t_i` in ``h`` and :math:`R_i^n` in ``E``: .. code-block:: python def convergence_rates(h, E): from math import log r = [log(E[i]/E[i-1])/log(h[i]/h[i-1]) for i in range(1, len(h))] return r Calling ``r_R_I = convergence_rates(dt, R_I)`` computes the sequence of rates :math:`r_0,r_1,\ldots,r_{m-2}` for the model :math:`R_I\sim\Delta t^r`, while the statements .. code-block:: python R = np.array(R) # two-dim. numpy array r_R = [convergence_rates(dt, R[:,n])[-1] for n in range(len(t_coarse))] compute the final rate :math:`r_{m-2}` for :math:`R^n\sim\Delta t^r` at each mesh point :math:`t_n` in the coarsest mesh. This latter computation deserves more explanation. Since ``R[i][n]`` holds the estimated truncation error :math:`R_i^n` on mesh :math:`i`, at point :math:`t_n` in the coarsest mesh, ``R[:,n]`` picks out the sequence :math:`R_i^n` for :math:`i=0,\ldots,m-1`. The ``convergence_rate`` function computes the rates at :math:`t_n`, and by indexing ``[-1]`` on the returned array from ``convergence_rate``, we pick the rate :math:`r_{m-2}`, which we believe is the best estimation since it is based on the two finest meshes. The ``estimate`` function is available in a module `trunc_empir.py `__. Let us apply this function to estimate the truncation error of the Forward Euler scheme. We need a function ``decay_FE(dt, N)`` that can compute :ref:`(688) ` at the points in a mesh with time step ``dt`` and ``N`` intervals: .. code-block:: python import numpy as np import trunc_empir def decay_FE(dt, N): dt = float(dt) t = np.linspace(0, N*dt, N+1) u_e = I*np.exp(-a*t) # exact solution, I and a are global u = u_e # naming convention when writing up the scheme R = np.zeros(N) for n in range(0, N): R[n] = (u[n+1] - u[n])/dt + a*u[n] # Theoretical expression for the trunction error R_a = 0.5*I*(-a)**2*np.exp(-a*t)*dt return R, t[:-1], R_a[:-1] if __name__ == '__main__': I = 1; a = 2 # global variables needed in decay_FE trunc_empir.estimate(decay_FE, T=2.5, N_0=6, m=4, makeplot=True) The estimated rates for the integrated truncation error :math:`R_I` become 1.1, 1.0, and 1.0 for this sequence of four meshes. All the rates for :math:`R^n`, computed as ``r_R``, are also very close to 1 at all mesh points. The agreement between the theoretical formula :ref:`(683) ` and the computed quantity (ref:ref:`(688) `) is very good, as illustrated in Figures :ref:`trunc:fig:FE:rates` and :ref:`trunc:fig:FE:error`. The program `trunc_decay_FE.py `__ was used to perform the simulations and it can easily be modified to test other schemes (see also :ref:`trunc:exer:decay:estimate`). .. _trunc:fig:FE:rates: .. figure:: R_series.png :width: 400 *Estimated truncation error at mesh points for different meshes* .. _trunc:fig:FE:error: .. figure:: R_error.png :width: 400 *Difference between theoretical and estimated truncation error at mesh points for different meshes* .. _trunc:decay:corr: Increasing the accuracy by adding correction terms -------------------------------------------------- .. index:: correction terms .. index:: single: truncation error; correction terms Now we ask the question: can we add terms in the differential equation that can help increase the order of the truncation error? To be precise, let us revisit the Forward Euler scheme for :math:`u'=-au`, insert the exact solution :math:`{u_{\small\mbox{e}}}`, include a residual :math:`R`, but also include new terms :math:`C`: .. _Eq:trunc:decay:FE:corr: .. math:: \tag{692} \lbrack D_t^+ {u_{\small\mbox{e}}} + a{u_{\small\mbox{e}}} = C + R \rbrack^n{\thinspace .} Inserting the Taylor expansions for :math:`[D_t^+{u_{\small\mbox{e}}}]^n` and keeping terms up to 3rd order in :math:`\Delta t` gives the equation .. math:: \frac{1}{2}{u_{\small\mbox{e}}}''(t_n)\Delta t - \frac{1}{6}{u_{\small\mbox{e}}}'''(t_n)\Delta t^2 + \frac{1}{24}{u_{\small\mbox{e}}}''''(t_n)\Delta t^3 + {\mathcal{O}(\Delta t^4)} = C^n + R^n{\thinspace .} Can we find :math:`C^n` such that :math:`R^n` is :math:`{\mathcal{O}(\Delta t^2)}`? Yes, by setting .. math:: C^n = \frac{1}{2}{u_{\small\mbox{e}}}''(t_n)\Delta t, we manage to cancel the first-order term and .. math:: R^n = \frac{1}{6}{u_{\small\mbox{e}}}'''(t_n)\Delta t^2 + {\mathcal{O}(\Delta t^3)}{\thinspace .} The correction term :math:`C^n` introduces :math:`\frac{1}{2}\Delta t u''` in the discrete equation, and we have to get rid of the derivative :math:`u''`. One idea is to approximate :math:`u''` by a second-order accurate finite difference formula, :math:`u''\approx (u^{n+1}-2u^n+u^{n-1})/\Delta t^2`, but this introduces an additional time level with :math:`u^{n-1}`. Another approach is to rewrite :math:`u''` in terms of :math:`u'` or :math:`u` using the ODE: .. math:: u'=-au\quad\Rightarrow\quad u''=-au' = -a(-au)= a^2u{\thinspace .} This means that we can simply set :math:`C^n = {\frac{1}{2}}a^2\Delta t u^n`. We can then either solve the discrete equation .. _Eq:trunc:decay:corr:FE:discrete: .. math:: \tag{693} [D_t^+ u = -au + {\frac{1}{2}}a^2\Delta t u]^n, or we can equivalently discretize the perturbed ODE .. _Eq:trunc:decay:corr:FE:ODE: .. math:: \tag{694} u' = -\hat au ,\quad \hat a = a(1 - {\frac{1}{2}}a\Delta t), by a Forward Euler method. That is, we replace the original coefficient :math:`a` by the perturbed coefficient :math:`\hat a`. Observe that :math:`\hat a\rightarrow a` as :math:`\Delta t\rightarrow 0`. The Forward Euler method applied to :ref:`(694) ` results in .. math:: [D_t^+ u = -a(1 - {\frac{1}{2}}a\Delta t)u]^n{\thinspace .} We can control our computations and verify that the truncation error of the scheme above is indeed :math:`{\mathcal{O}(\Delta t^2)}`. Another way of revealing the fact that the perturbed ODE leads to a more accurate solution is to look at the amplification factor. Our scheme can be written as .. math:: u^{n+1} = Au^n,\quad A = 1-\hat a\Delta t = 1 - p + {\frac{1}{2}}p^2,\quad p=a\Delta t, The amplification factor :math:`A` as a function of :math:`p=a\Delta t` is seen to be the first three terms of the Taylor series for the exact amplification factor :math:`e^{-p}`. The Forward Euler scheme for :math:`u=-au` gives only the first two terms :math:`1-p` of the Taylor series for :math:`e^{-p}`. That is, using :math:`\hat a` increases the order of the accuracy in the amplification factor. Instead of replacing :math:`u''` by :math:`a^2u`, we use the relation :math:`u''=-au'` and add a term :math:`-{\frac{1}{2}}a\Delta t u'` in the ODE: .. math:: u'=-au - {\frac{1}{2}}a\Delta t u'\quad\Rightarrow\quad \left( 1 + {\frac{1}{2}}a\Delta t\right) u' = -au{\thinspace .} Using a Forward Euler method results in .. math:: \left( 1 + {\frac{1}{2}}a\Delta t\right)\frac{u^{n+1}-u^n}{\Delta t} = -au^n, which after some algebra can be written as .. math:: u^{n+1} = \frac{1 - {\frac{1}{2}}a\Delta t}{1+{\frac{1}{2}}a\Delta t}u^n{\thinspace .} This is the same formula as the one arising from a Crank-Nicolson scheme applied to :math:`u'=-au`! It now recommended to do :ref:`trunc:exer:decay:corr:BE` and repeat the above steps to see what kind of correction term is needed in the Backward Euler scheme to make it second order. The Crank-Nicolson scheme is a bit more challenging to analyze, but the ideas and techniques are the same. The discrete equation reads .. math:: [D_t u = -au ]^{n+\frac{1}{2}}, and the truncation error is defined through .. math:: [D_t {u_{\small\mbox{e}}} + a\overline{{u_{\small\mbox{e}}}}^{t} = C + R]^{n+\frac{1}{2}}, where we have added a correction term. We need to Taylor expand both the discrete derivative and the arithmetic mean with aid of :ref:`(658) `-:ref:`(659) ` and :ref:`(674) `-:ref:`(675) `, respectively. The result is .. math:: \frac{1}{24}{u_{\small\mbox{e}}}'''(t_{n+\frac{1}{2}})\Delta t^2 + {\mathcal{O}(\Delta t^4)} + \frac{a}{8}{u_{\small\mbox{e}}}''(t_{n+\frac{1}{2}})\Delta t^2 + {\mathcal{O}(\Delta t^4)} = C^{n+\frac{1}{2}} + R^{n+\frac{1}{2}}{\thinspace .} The goal now is to make :math:`C^{n+\frac{1}{2}}` cancel the :math:`\Delta t^2` terms: .. math:: C^{n+\frac{1}{2}} = \frac{1}{24}{u_{\small\mbox{e}}}'''(t_{n+\frac{1}{2}})\Delta t^2 + \frac{a}{8}{u_{\small\mbox{e}}}''(t_{n})\Delta t^2{\thinspace .} Using :math:`u'=-au`, we have that :math:`u''=a^2u`, and we find that :math:`u'''=-a^3u`. We can therefore solve the perturbed ODE problem .. math:: u' = -\hat a u,\quad \hat a = a(1 - \frac{1}{12}a^2\Delta t^2), by the Crank-Nicolson scheme and obtain a method that is of fourth order in :math:`\Delta t`. :ref:`trunc:exer:decay:corr:verify` encourages you to implement these correction terms and calculate empirical convergence rates to verify that higher-order accuracy is indeed obtained in real computations. Extension to variable coefficients (1) ----------------------------------------------- Let us address the decay ODE with variable coefficients, .. math:: u'(t) = -a(t)u(t) + b(t), discretized by the Forward Euler scheme, .. _Eq:_auto259: .. math:: \tag{695} [D_t^+ u = -au + b]^n {\thinspace .} The truncation error :math:`R` is as always found by inserting the exact solution :math:`{u_{\small\mbox{e}}}(t)` in the discrete scheme: .. _Eq:_auto260: .. math:: \tag{696} [D_t^+ {u_{\small\mbox{e}}} + a{u_{\small\mbox{e}}} - b = R]^n {\thinspace .} Using :ref:`(664) `-:ref:`(665) `, .. math:: {u_{\small\mbox{e}}}'(t_n) - \frac{1}{2}{u_{\small\mbox{e}}}''(t_n)\Delta t + {\mathcal{O}(\Delta t^2)} + a(t_n){u_{\small\mbox{e}}}(t_n) - b(t_n) = R^n {\thinspace .} Because of the ODE, .. math:: {u_{\small\mbox{e}}}'(t_n) + a(t_n){u_{\small\mbox{e}}}(t_n) - b(t_n) =0, so we are left with the result .. _Eq:trunc:decay:vc:R: .. math:: \tag{697} R^n = -\frac{1}{2}{u_{\small\mbox{e}}}''(t_n)\Delta t + {\mathcal{O}(\Delta t^2)} \ {\thinspace .} We see that the variable coefficients do not pose any additional difficulties in this case. :ref:`trunc:exer:decay:varcoeff:CN` takes the analysis above one step further to the Crank-Nicolson scheme. Exact solutions of the finite difference equations -------------------------------------------------- .. index:: verification Having a mathematical expression for the numerical solution is very valuable in program verification since we then know the exact numbers that the program should produce. Looking at the various formulas for the truncation errors in :ref:`(658) `-:ref:`(659) ` and :ref:`(678) `-:ref:`(679) ` in the section :ref:`trunc:table`, we see that all but two of the :math:`R` expressions contain a second or higher order derivative of :math:`{u_{\small\mbox{e}}}`. The exceptions are the geometric and harmonic means where the truncation error involves :math:`{u_{\small\mbox{e}}}'` and even :math:`{u_{\small\mbox{e}}}` in case of the harmonic mean. So, apart from these two means, choosing :math:`{u_{\small\mbox{e}}}` to be a linear function of :math:`t`, :math:`{u_{\small\mbox{e}}} = ct+d` for constants :math:`c` and :math:`d`, will make the truncation error vanish since :math:`{u_{\small\mbox{e}}}''=0`. Consequently, the truncation error of a finite difference scheme will be zero since the various approximations used will all be exact. This means that the linear solution is an exact solution of the discrete equations. In a particular differential equation problem, the reasoning above can be used to determine if we expect a linear :math:`{u_{\small\mbox{e}}}` to fulfill the discrete equations. To actually prove that this is true, we can either compute the truncation error and see that it vanishes, or we can simply insert :math:`{u_{\small\mbox{e}}}(t)=ct+d` in the scheme and see that it fulfills the equations. The latter method is usually the simplest. It will often be necessary to add some source term to the ODE in order to allow a linear solution. Many ODEs are discretized by centered differences. From the section :ref:`trunc:table` we see that all the centered difference formulas have truncation errors involving :math:`{u_{\small\mbox{e}}}'''` or higher-order derivatives. A quadratic solution, e.g., :math:`{u_{\small\mbox{e}}}(t) =t^2 + ct + d`, will then make the truncation errors vanish. This observation can be used to test if a quadratic solution will fulfill the discrete equations. Note that a quadratic solution will not obey the equations for a Crank-Nicolson scheme for :math:`u'=-au+b` because the approximation applies an arithmetic mean, which involves a truncation error with :math:`{u_{\small\mbox{e}}}''`. .. _trunc:decay:gen: Computing truncation errors in nonlinear problems ------------------------------------------------- The general nonlinear ODE .. _Eq:trunc:decay:gen:ode: .. math:: \tag{698} u'=f(u,t), can be solved by a Crank-Nicolson scheme .. _Eq:trunc:decay:gen:ode:fdm: .. math:: \tag{699} [D_t u=\overline{f}^{t}]^{n+\frac{1}{2}}{\thinspace .} The truncation error is as always defined as the residual arising when inserting the exact solution :math:`{u_{\small\mbox{e}}}` in the scheme: .. _Eq:trunc:decay:gen:ode:CN: .. math:: \tag{700} [D_t {u_{\small\mbox{e}}} - \overline{f}^{t}= R]^{n+\frac{1}{2}}{\thinspace .} Using :ref:`(674) `-:ref:`(675) ` for :math:`\overline{f}^{t}` results in .. math:: \begin{align*} [\overline{f}^{t}]^{n+\frac{1}{2}} &= \frac{1}{2}(f({u_{\small\mbox{e}}}^n,t_n) + f({u_{\small\mbox{e}}}^{n+1},t_{n+1}))\\ &= f({u_{\small\mbox{e}}}^{n+\frac{1}{2}},t_{n+\frac{1}{2}}) + \frac{1}{8}{u_{\small\mbox{e}}}''(t_{n+\frac{1}{2}})\Delta t^2 + {\mathcal{O}(\Delta t^4)}{\thinspace .} \end{align*} With :ref:`(658) `-:ref:`(659) ` the discrete equations :ref:`(700) ` lead to .. math:: {u_{\small\mbox{e}}}'(t_{n+\frac{1}{2}}) + \frac{1}{24}{u_{\small\mbox{e}}}'''(t_{n+\frac{1}{2}})\Delta t^2 - f({u_{\small\mbox{e}}}^{n+\frac{1}{2}},t_{n+\frac{1}{2}}) - \frac{1}{8}{u_{\small\mbox{e}}}''(t_{n+\frac{1}{2}})\Delta t^2 + {\mathcal{O}(\Delta t^4)} = R^{n+\frac{1}{2}}{\thinspace .} Since :math:`{u_{\small\mbox{e}}}'(t_{n+\frac{1}{2}}) - f({u_{\small\mbox{e}}}^{n+\frac{1}{2}},t_{n+\frac{1}{2}})=0`, the truncation error becomes .. math:: R^{n+\frac{1}{2}} = (\frac{1}{24}{u_{\small\mbox{e}}}'''(t_{n+\frac{1}{2}}) - \frac{1}{8}{u_{\small\mbox{e}}}''(t_{n+\frac{1}{2}})) \Delta t^2{\thinspace .} The computational techniques worked well even for this nonlinear ODE.