.. !split

.. _ch:formulas:

Appendix: Useful formulas
%%%%%%%%%%%%%%%%%%%%%%%%%

.. _sec:form:fdop:

Finite difference operator notation
===================================

.. _Eq:Dop:fd1:center:

.. math::

    \tag{609}
    u'(t_n) \approx
        \lbrack D_tu\rbrack^n = \frac{u^{n+\frac{1}{2}} - u^{n-\frac{1}{2}}}{\Delta t} 
        

.. _Eq:Dop:fd1:center2:

.. math::

    \tag{610}
    u'(t_n) \approx
        \lbrack D_{2t}u\rbrack^n = \frac{u^{n+1} - u^{n-1}}{2\Delta t}
        
        

.. _Eq:Dop:fd1:bw:

.. math::

    \tag{611}
    u'(t_n) =
        \lbrack D_t^-u\rbrack^n = \frac{u^{n} - u^{n-1}}{\Delta t}
        
        

.. _Eq:Dop:fd1:fw:

.. math::

    \tag{612}
    u'(t_n) \approx
        \lbrack D_t^+u\rbrack^n = \frac{u^{n+1} - u^{n}}{\Delta t}
        
        

.. _Eq:Dop:fd1:theta:

.. math::

    \tag{613}
    u'(t_{n+\theta}) =
        \lbrack \bar D_tu\rbrack^{n+\theta} = \frac{u^{n+1} - u^{n}}{\Delta t}
        
        

.. _Eq:Dop:fd1:bw2:

.. math::

    \tag{614}
    u'(t_n) \approx
        \lbrack D_t^{2-}u\rbrack^n = \frac{3u^{n} - 4u^{n-1} + u^{n-2}}{2\Delta t}
        
        

.. _Eq:Dop:fd2:center:

.. math::

    \tag{615}
    u''(t_n) \approx
        \lbrack D_tD_t u\rbrack^n = \frac{u^{n+1} - 2u^{n} + u^{n-1}}{\Delta t^2}
        
        

.. _Eq:Dop:avg:arith:

.. math::

    \tag{616}
    u(t_{n+\frac{1}{2}}) \approx \lbrack \overline{u}^{t}\rbrack^{n+\frac{1}{2}}
        = \frac{1}{2}(u^{n+1} + u^n)
        
        

.. _Eq:Dop:avg:geom:

.. math::

    \tag{617}
    u(t_{n+\frac{1}{2}})^2 \approx \lbrack \overline{u^2}^{t,g}\rbrack^{n+\frac{1}{2}}
        = u^{n+1}u^n
        
        

.. _Eq:Dop:avg:harm:

.. math::

    \tag{618}
    u(t_{n+\frac{1}{2}}) \approx \lbrack \overline{u}^{t,h}\rbrack^{n+\frac{1}{2}}
        = \frac{2}{\frac{1}{u^{n+1}} + \frac{1}{u^n}}
        
        

.. _Eq:Dop:avg:theta:

.. math::

    \tag{619}
    u(t_{n+\theta}) \approx \lbrack \overline{u}^{t,\theta}\rbrack^{n+\theta}
        = \theta u^{n+1} + (1-\theta)u^n ,
        
        

.. _Eq:_auto253:

.. math::

    \tag{620}
    \qquad t_{n+\theta}=\theta t_{n+1} + (1-\theta)t_{n-1}
        
        

Some may wonder why :math:`\theta` is absent on the right-hand side
of :ref:`(613) <Eq:Dop:fd1:theta>`. The fraction is an approximation to the
derivative at the point :math:`t_{n+\theta}=\theta t_{n+1} + (1-theta) t_{n}`.

.. _sec:form:truncerr:

Truncation errors of finite difference approximations
=====================================================

.. math::
        
        {u_{\small\mbox{e}}}'(t_n) =
        [D_t{u_{\small\mbox{e}}}]^n + R^n = \frac{{u_{\small\mbox{e}}}^{n+\frac{1}{2}} - {u_{\small\mbox{e}}}^{n-\frac{1}{2}}}{\Delta t} +R^n\nonumber,
        

.. _Eq:form:trunc:fd1:center:

.. math::

    \tag{621}
    R^n = -\frac{1}{24}{u_{\small\mbox{e}}}'''(t_n)\Delta t^2 + {\cal O}(\Delta t^4)
        
        

.. math::
          
        {u_{\small\mbox{e}}}'(t_n) =
        [D_{2t}{u_{\small\mbox{e}}}]^n +R^n = \frac{{u_{\small\mbox{e}}}^{n+1} - {u_{\small\mbox{e}}}^{n-1}}{2\Delta t} +
        R^n\nonumber,
        

.. _Eq:form:trunc:fd1:center2:

.. math::

    \tag{622}
    R^n = -\frac{1}{6}{u_{\small\mbox{e}}}'''(t_n)\Delta t^2 + {\cal O}(\Delta t^4)
        
        

.. math::
          
        {u_{\small\mbox{e}}}'(t_n) =
        [D_t^-{u_{\small\mbox{e}}}]^n +R^n = \frac{{u_{\small\mbox{e}}}^{n} - {u_{\small\mbox{e}}}^{n-1}}{\Delta t}
        +R^n\nonumber,
        

.. _Eq:form:trunc:fd1:bw:

.. math::

    \tag{623}
    R^n = -\frac{1}{2}{u_{\small\mbox{e}}}''(t_n)\Delta t + {\cal O}(\Delta t^2)
        
        

.. math::
          
        {u_{\small\mbox{e}}}'(t_n) =
        [D_t^+{u_{\small\mbox{e}}}]^n +R^n = \frac{{u_{\small\mbox{e}}}^{n+1} - {u_{\small\mbox{e}}}^{n}}{\Delta t}
        +R^n\nonumber,
        

.. _Eq:form:trunc:fd1:fw:

.. math::

    \tag{624}
    R^n = \frac{1}{2}{u_{\small\mbox{e}}}''(t_n)\Delta t + {\cal O}(\Delta t^2)
        
        

.. math::
          
        {u_{\small\mbox{e}}}'(t_{n+\theta}) =
        [\bar D_t{u_{\small\mbox{e}}}]^{n+\theta} +R^{n+\theta} = \frac{{u_{\small\mbox{e}}}^{n+1} - {u_{\small\mbox{e}}}^{n}}{\Delta t}
        +R^{n+\theta}\nonumber,
        

.. math::
          
        R^{n+\theta} = -\frac{1}{2}(1-2\theta){u_{\small\mbox{e}}}''(t_{n+\theta})\Delta t +
        \frac{1}{6}((1 - \theta)^3 - \theta^3){u_{\small\mbox{e}}}'''(t_{n+\theta})\Delta t^2 +
        \nonumber
        

.. _Eq:form:trunc:fd1:theta:

.. math::

    \tag{625}
    \quad {\cal O}(\Delta t^3)
        
        

.. math::
          
        {u_{\small\mbox{e}}}'(t_n) =
        [D_t^{2-}{u_{\small\mbox{e}}}]^n +R^n = \frac{3{u_{\small\mbox{e}}}^{n} - 4{u_{\small\mbox{e}}}^{n-1} + {u_{\small\mbox{e}}}^{n-2}}{2\Delta t}
        +R^n\nonumber,
        

.. _Eq:form:trunc:fd1:bw2:

.. math::

    \tag{626}
    R^n = \frac{1}{3}{u_{\small\mbox{e}}}'''(t_n)\Delta t^2 + {\cal O}(\Delta t^3)
        
        

.. math::
          
        {u_{\small\mbox{e}}}''(t_n) =
        [D_tD_t {u_{\small\mbox{e}}}]^n +R^n = \frac{{u_{\small\mbox{e}}}^{n+1} - 2{u_{\small\mbox{e}}}^{n} + {u_{\small\mbox{e}}}^{n-1}}{\Delta t^2}
        +R^n\nonumber,
        

.. _Eq:form:trunc:fd2:center:

.. math::

    \tag{627}
    R^n = -\frac{1}{12}{u_{\small\mbox{e}}}''''(t_n)\Delta t^2 + {\cal O}(\Delta t^4)
        
        

.. math::
        
        {u_{\small\mbox{e}}}(t_{n+\theta}) = [\overline{{u_{\small\mbox{e}}}}^{t,\theta}]^{n+\theta} +R^{n+\theta}
        = \theta {u_{\small\mbox{e}}}^{n+1} + (1-\theta){u_{\small\mbox{e}}}^n +R^{n+\theta},\nonumber
        

.. _Eq:form:trunc:avg:theta:

.. math::

    \tag{628}
    R^{n+\theta} = -\frac{1}{2}{u_{\small\mbox{e}}}''(t_{n+\theta})\Delta t^2\theta (1-\theta) +
        {\cal O}(\Delta t^3)
        {\thinspace .}
        
        

.. _sec:form:fdexp:

Finite differences of exponential functions
===========================================

Complex exponentials
~~~~~~~~~~~~~~~~~~~~

Let :math:`u^n = \exp{(i\omega n\Delta t)} = e^{i\omega t_n}`.

.. _Eq:form:exp:fd2:center:

.. math::

    \tag{629}
    [D_tD_t u]^n = u^n \frac{2}{\Delta t}(\cos \omega\Delta t - 1) =
        -\frac{4}{\Delta t}\sin^2\left(\frac{\omega\Delta t}{2}\right),
        
        

.. _Eq:form:exp:fd1:fw:

.. math::

    \tag{630}
    [D_t^+ u]^n = u^n\frac{1}{\Delta t}(\exp{(i\omega\Delta t)} - 1),
        
        

.. _Eq:form:exp:fd1:bw:

.. math::

    \tag{631}
    [D_t^- u]^n = u^n\frac{1}{\Delta t}(1 - \exp{(-i\omega\Delta t)}),
        
        

.. _Eq:form:exp:fd1:center:

.. math::

    \tag{632}
    [D_t u]^n = u^n\frac{2}{\Delta t}i\sin{\left(\frac{\omega\Delta t}{2}\right)},
        
        

.. _Eq:form:exp:fd1c:center:

.. math::

    \tag{633}
    [D_{2t} u]^n = u^n\frac{1}{\Delta t}i\sin{(\omega\Delta t)}
        
        {\thinspace .}
        

Real exponentials
~~~~~~~~~~~~~~~~~

Let :math:`u^n = \exp{(\omega n\Delta t)} = e^{\omega t_n}`.

.. _Eq:form:rexp:fd2:center:

.. math::

    \tag{634}
    [D_tD_t u]^n = u^n \frac{2}{\Delta t}(\cos \omega\Delta t - 1) =
        -\frac{4}{\Delta t}\sin^2\left(\frac{\omega\Delta t}{2}\right),
        
        

.. _Eq:form:rexp:fd1:fw:

.. math::

    \tag{635}
    [D_t^+ u]^n = u^n\frac{1}{\Delta t}(\exp{(i\omega\Delta t)} - 1),
        
        

.. _Eq:form:rexp:fd1:bw:

.. math::

    \tag{636}
    [D_t^- u]^n = u^n\frac{1}{\Delta t}(1 - \exp{(-i\omega\Delta t)}),
        
        

.. _Eq:form:rexp:fd1:center:

.. math::

    \tag{637}
    [D_t u]^n = u^n\frac{2}{\Delta t}i\sin{\left(\frac{\omega\Delta t}{2}\right)},
        
        

.. _Eq:form:rexp:fd1c:center:

.. math::

    \tag{638}
    [D_{2t} u]^n = u^n\frac{1}{\Delta t}i\sin{(\omega\Delta t)}
        
        {\thinspace .}
        

.. _sec:form:fdtn:

Finite differences of :math:`t^n`
=================================

The following results are useful when checking if a polynomial term in a
solution fulfills the discrete equation for the numerical method.

.. _Eq:form:Dop:tn:fw:

.. math::

    \tag{639}
    \lbrack D_t^+ t\rbrack^n = 1,
        
        

.. _Eq:form:Dop:tn:bw:

.. math::

    \tag{640}
    \lbrack D_t^- t\rbrack^n = 1,
        
        

.. _Eq:form:Dop:tn:cn:

.. math::

    \tag{641}
    \lbrack D_t t\rbrack^n = 1,
        
        

.. _Eq:form:Dop:tn:2cn:

.. math::

    \tag{642}
    \lbrack D_{2t} t\rbrack^n = 1,
        
        

.. _Eq:form:Dop2:tn:cn:

.. math::

    \tag{643}
    \lbrack D_{t}D_t t\rbrack^n = 0
        
        {\thinspace .}
        

The next formulas concern the action of difference operators on a :math:`t^2` term.

.. _Eq:form:Dop:tn2:fw:

.. math::

    \tag{644}
    \lbrack D_t^+ t^2\rbrack^n = (2n+1)\Delta t,
        
        

.. _Eq:form:Dop:tn2:bw:

.. math::

    \tag{645}
    \lbrack D_t^- t^2\rbrack^n = (2n-1)\Delta t,
        
        

.. _Eq:form:Dop:tn2:cn:

.. math::

    \tag{646}
    \lbrack D_t t^2\rbrack^n = 2n\Delta t,
        
        

.. _Eq:form:Dop:tn2:2cn:

.. math::

    \tag{647}
    \lbrack D_{2t} t^2\rbrack^n = 2n\Delta t,
        
        

.. _Eq:form:Dop2:tn2:cn:

.. math::

    \tag{648}
    \lbrack D_{t}D_t t^2\rbrack^n = 2,
        
        

Finally, we present formulas for a :math:`t^3` term:

.. _Eq:form:Dop:tn3:fw:

.. math::

    \tag{649}
    \lbrack D_t^+ t^3\rbrack^n = 3(n\Delta t)^2 + 3n\Delta t^2 + \Delta t^2,
        
        

.. _Eq:form:Dop:tn3:bw:

.. math::

    \tag{650}
    \lbrack D_t^- t^3\rbrack^n = 3(n\Delta t)^2 - 3n\Delta t^2 + \Delta t^2,
        
        

.. _Eq:form:Dop:tn3:cn:

.. math::

    \tag{651}
    \lbrack D_t t^3\rbrack^n = 3(n\Delta t)^2 + \frac{1}{4}\Delta t^2,
        
        

.. _Eq:form:Dop:tn3:2cn:

.. math::

    \tag{652}
    \lbrack D_{2t} t^3\rbrack^n = 3(n\Delta t)^2 + \Delta t^2,
        
        

.. _Eq:form:Dop2:tn3:cn:

.. math::

    \tag{653}
    \lbrack D_{t}D_t t^3\rbrack^n = 6n\Delta t,
        
        

Software
--------

Application of finite difference operators to polynomials and exponential
functions, resulting in the formulas above, can easily be computed by
some ``sympy`` code (from the file `lib.py <http://tinyurl.com/nu656p2/formulas/lib.py>`__):

.. code-block:: python

    from sympy import *
    t, dt, n, w = symbols('t dt n w', real=True)
    
    # Finite difference operators
    
    def D_t_forward(u):
        return (u(t + dt) - u(t))/dt
    
    def D_t_backward(u):
        return (u(t) - u(t-dt))/dt
    
    def D_t_centered(u):
        return (u(t + dt/2) - u(t-dt/2))/dt
    
    def D_2t_centered(u):
        return (u(t + dt) - u(t-dt))/(2*dt)
    
    def D_t_D_t(u):
        return (u(t + dt) - 2*u(t) + u(t-dt))/(dt**2)
    
    
    op_list = [D_t_forward, D_t_backward,
               D_t_centered, D_2t_centered, D_t_D_t]
    
    def ft1(t):
        return t
    
    def ft2(t):
        return t**2
    
    def ft3(t):
        return t**3
    
    def f_expiwt(t):
        return exp(I*w*t)
    
    def f_expwt(t):
        return exp(w*t)
    
    func_list = [ft1, ft2, ft3, f_expiwt, f_expwt]

To see the results, one can now make a simple loop over the different
type of functions and the various operators associated with them:

.. code-block:: python

    for func in func_list:
        for op in op_list:
            f = func
            e = op(f)
            e = simplify(expand(e))
            print e
            if func in [f_expiwt, f_expwt]:
                e = e/f(t)
            e = e.subs(t, n*dt)
            print expand(e)
            print factor(simplify(expand(e)))