Computation and visualization of truncation error and round-off error
​by Sopanant Datta and Taweetham Limpanuparb

Introduction

Round-off error and truncation error are common inherent sources of error in numerical computation. This post shows the computation and visualization of round-off and truncation errors in specific examples.

Round-off error

A round-off error is the difference between a calculated (approximated with finite-precision arithmetic) value and a reference (exact) value. One area whereby this source of error can be a crucial consideration is the computation of Associated Legendre polynomials (ALP) and spherical harmonics (SH) in quantum chemistry applications. The article “Associated Legendre Polynomials and Spherical Harmonics Computation for Chemistry Applications” (Limpanuparb & Milthorpe, 2014), with the objective of developing an algorithm to efficiently and accurately compute these functions, presents the computation and visualization of round-off errors in ALPs and SHs.

Truncation error

A truncation error is the difference between a numerical approximation of a value and the analytical or true value. An example of a truncation error arises from the approximation of an integral by a sum (truncated infinite series). The article “A Remarkable Identity Involving Bessel Functions” (Dominici, Gill & Limpanuparb, 2012) presents and proves an identity of Bessel function integrals and sums and demonstrates the computation and visualization of truncation errors in the identity.

Computation and visualization

Round-off error

Definitions of calculated values, reference values and errors

◼
  • Calculated values
    X
    cal
    are computed using an algorithm developed by the authors described in the article. These values are imported for further calculations of errors.
    • ◼
  • Reference values
    X
    ref
    of ALPs are computed using the function
    P
    l_,m_
    [x_] , where x = cosθ, and degree l and order m satisfy the conditions l ≥ 0 and -l ≤ m ≤ l.
    • P
    l_,m_
    [x_]:=
    (2l+1)(l-m)!
    2π(l+m)!
    LegendreP[l,m,x]
    ◼
  • Absolute error
    ϵ
    a
    andrelativeerror
    ϵ
    r
    are defined as follows:​
    ϵ
    a
    =
    max(
    log
    10
    |
    X
    cal
    -
    X
    ref
    |,-324
    )​
    ϵ
    r
    =
    minmax
    log
    10
    
    X
    cal
    X
    ref
    -1,-16
    ),0).

    • Contour plot of magnitude of errors​Example of evaluation for l ≤ 100 and θ =
    π
    100

    A list of
    X
    cal
    values is imported (cval1) and
    X
    ref
    values are computed (rval1). The magnitudes of
    X
    cal
    and
    X
    ref
    values are calculated and contained in lists mcval1 and mrval1, respectively.
    θ=
    π
    100
    ;​​cval1=Import[NotebookDirectory[]<>"Code/cPlm_100_PI100_0.out","Table"];(*
    X
    cal
    forθ=0*)Length[cval1]-1(*=100*)​​rval1=Table[N[
    P
    l,m
    [Cos[θ]],100],{l,0,100},{m,0,l}];LL=Length[rval1]-1(*=100*)​​mcval1=Table[If[m≤l,Log10[Abs[cval1〚l+1〛〚m+1〛]]/.{Indeterminate-323}],{l,0,LL},{m,0,LL}];​​mrval1=Table[If[m≤l,Log10[Abs[rval1〚l+1〛〚m+1〛]]/.{-∞-323}],{l,0,LL},{m,0,LL}];
    Out[]=
    100
    Out[]=
    100
    The magnitudes of
    ϵ
    a
    and
    ϵ
    r
    are computed and contained in lists aepad and repad, respectively.
    In[]:=
    ae=Log10[Abs[cval1-rval1]]/.{Indeterminate-323};​​ae1=Table[If[m≤l,ae〚l+1〛〚m+1〛],{l,0,LL},{m,0,LL}];​​re1=Table[Which[m≤l&&rval1〚l+1〛〚m+1〛≠0,Max[Log10[Abs[cval1〚l+1〛〚m+1〛/rval1〚l+1〛〚m+1〛-1]],-16],m≤l,-16,True,],{l,0,LL},{m,0,LL}];
    The magnitudes of
    X
    cal
    ,
    X
    ref
    ,
    ϵ
    a
    and
    ϵ
    r
    can be visualized with a contour plot. ListContourPlot is used as an example here. The plots of
    X
    cal
    and
    X
    ref
    magnitudes appear to be identical. This is supported by how
    ϵ
    a
    and
    ϵ
    r
    values are all under
    -10
    10
    , as shown in the plots of
    ϵ
    a
    and
    ϵ
    r
    magnitudes.
    In[]:=
    mcplot=ListContourPlot[mcval1,FrameLabel{m,l},PlotLegendsAutomatic,DataRange{{0,100},{0,100}},PlotLabel->"Magnitude of
    X
    cal
    "];​​mrplot=ListContourPlot[mrval1,FrameLabel{m,l},PlotLegendsAutomatic,DataRange{{0,100},{0,100}},PlotLabel->"Magnitude of
    X
    ref
    "];​​aeplot=ListContourPlot[ae1,FrameLabel{m,l},PlotLegendsAutomatic,DataRange{{0,100},{0,100}},PlotLabel->"Magnitude of
    ϵ
    a
    "];​​replot=ListContourPlot[re1,FrameLabel{m,l},PlotLegendsAutomatic,DataRange{{0,100},{0,100}},PlotLabel->"Magnitude of
    ϵ
    r
    "];​​GraphicsRow[{mcplot,mrplot,aeplot,replot},ImageSize1200]
    Out[]=

    Scatter plot of magnitude of errors​Example of evaluation for l ≤ 1000 and θ =
    0,
    π
    100
    ,
    π
    4
    ,
    49π
    100
    ,
    π
    2

    A list of
    X
    cal
    values is imported (cval2) and
    X
    ref
    values are computed (rval2). The magnitudes of
    X
    cal
    ,
    X
    ref
    ,
    ϵ
    a
    and
    ϵ
    r
    values are calculated and contained in lists mcval2, mrval2, ae2 and re2, respectively.
    cval2=Table[Import[NotebookDirectory[]<>"Code/"<>file,"Table"],{file,{"cPlm_1000_0_0.out","cPlm_1000_PI100_0.out","cPlm_1000_PI4_0.out","cPlm_1000_49PI100_0.out","cPlm_1000_PI2_0.out"}}];​​rval2=TableN[
    P
    1000,m
    [Cos[θ]],1000],θ,0,
    π
    100
    ,
    π
    4
    ,
    49π
    100
    ,
    π
    2
    ,{l,1000,1000},{m,0,l};​​mcval2=Table[Log10[Abs[cval2〚i〛〚1〛]]/.{Indeterminate-323},{i,Length[cval2]}];​​mrval2=Table[Log10[Abs[rval2〚i〛〚1〛]]/.{Indeterminate-323},{i,Length[rval2]}];​​ae2=Table[Log10[Abs[cval2〚i〛〚1〛-rval2〚i〛〚1〛]]/.{Indeterminate-323},{i,Length[cval2]}];​​re2=Table[If[rval2〚i〛〚1〛〚m+1〛≠0,Max[Log10[Abs[cval2〚i〛〚1〛〚m+1〛/rval2〚i〛〚1〛〚m+1〛-1]],-16]],{i,Length[cval2]},{m,0,LL}];
    The magnitudes of
    X
    cal
    ,
    X
    ref
    ,
    ϵ
    a
    and
    ϵ
    r
    are visualized using ListPlot.
    In[]:=
    mcplot=ListPlotmcval2,PlotRange{0,-20},PlotLegends"θ=0","θ=
    π
    100
    ","θ=
    π
    4
    ","θ=
    49
    100
    π","θ=
    π
    2
    ",AxesLabel{m,"
    χ
    cal
    "},DataRange{0,1000};​​mrplot=ListPlotmrval2,PlotRange{0,-20},PlotLegends"θ=0","θ=
    π
    100
    ","θ=
    π
    4
    ","θ=
    49
    100
    π","θ=
    π
    2
    ",AxesLabel{m,"
    χ
    ref
    "},DataRange{0,1000};​​aeplot=ListPlotae2,PlotRange{0,-20},PlotLegends"θ=0","θ=
    π
    100
    ","θ=
    π
    4
    ","θ=
    49
    100
    π","θ=
    π
    2
    ",AxesLabel{m,"
    ϵ
    a
    "},DataRange{0,1000};​​replot=ListPlotre2,PlotRange{0,-20},PlotLegends"θ=0","θ=
    π
    100
    ","θ=
    π
    4
    ","θ=
    49
    100
    π","θ=
    π
    2
    ",AxesLabel{m,"
    ϵ
    r
    "},DataRange{0,1000};​​GraphicsRow[{mcplot,mrplot,aeplot,replot},Spacings->0,ImageSize->1300]

    Truncation error

    Integral and truncated sum

    The followings are the integral and truncated sum of Bessel functions, respectively.​
    ∞
    ∫
    0
    J
    5
    2
    (at)
    J
    5
    2
    (bt)
    t
    t
    =
    1
    5
    (|a|)
    a
    5
    2
    (|b|)
    b
    
    (min(|a|,|b|)
    (max(|a|,|b|)
    
    ​​
    T
    N
    (a,b)
    =
    N
    ∑
    n=1
    J
    5
    2
    (an)
    J
    5
    2
    (bn)
    n
    The integral (X; solid) and truncated sum
    T
    10
    (a,b)
    (Y; dashed) can be defined and plotted as follows .
    In[]:=
    J
    l_
    [x_]:=BesselJ[l,x];​​X=Table
    1
    5
    5/2
    Min[a,b]
    Max[a,b]
    ,a,
    π
    4
    ,
    3
    4
    π,
    π
    4
    ;​​Y=Table
    10
    ∑
    n=1
    J
    5
    2
    [an]
    J
    5
    2
    [bn]
    n
    ,a,
    π
    4
    ,
    3
    4
    π,
    π
    4
    ;
    In[]:=
    Z=Union[Y,X];​​PlotZ,{b,0,2π},PlotStyle{{Darker[Green],Thickness[0.002]},{Darker[Green],Thickness[0.002]},{Darker[Green],Thickness[0.002]},{RGBColor[1,0,0],Thickness[0.0015],Dashed},{RGBColor[1,0,0],Thickness[0.0015],Dashed},{RGBColor[1,0,0],Thickness[0.0015],Dashed},{RGBColor[0,0,1],Thickness[0.0015],DotDashed},{RGBColor[0,0,1],Thickness[0.0015],DotDashed},{RGBColor[0,0,1],Thickness[0.0015],DotDashed}},PlotRange{0,.2},AxesLabel{b,""},​​Prolog​​Text[Style["a",FontSize9,Italic],{1.6,0.04},{-1,0}],​​TextStyle"=
    π
    4
    ",FontSize9,{1.75,0.039},{-1,0},​​Text[Style["a",FontSize9,Italic],{2.5,.07},{-1,0}],​​TextStyle"=
    2π
    4
    ",FontSize9,{2.65,.07},{-1,0},​​Text[Style["a",FontSize9,Italic],{3.1,.12},{-1,0}],​​TextStyle"=
    3π
    4
    ",FontSize9,{3.25,.12},{-1,0}​​,​​Ticks
    π
    4
    ,
    2π
    4
    ,
    3π
    4
    ,
    4π
    4
    ,
    5π
    4
    ,
    6π
    4
    ,
    7π
    4
    ,
    8π
    4
    ,Automatic​​​​

    Visualization of error