5
$\begingroup$

Consider a data set with some numerical uncertainty in the y-values, such as the following

g[x_] := x^3
data = Table[{x, Around[g[x], RandomReal[{0, .1}]]}, {x, 0, 2, .1}]

I want to numerically integrate the data and find the uncertainty in the result. The only way I know how to numerically integrate the data in Mathematica is to first interpolate, such as

interp = Interpolation[data]
int[y_] := NIntegrate[interp[x], {x, 0, y}]

When I do so with the data with uncertainty, however, I get

The integrand InterpolatingFunction[...] has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,0.1}}

Is Mathematica capable of doing these numerical integrals with uncertainties? And if so, is it documented how it estimates the uncertainty? If not, is there a standard method/workaround to this (I guess something like Monte Carlo resampling?)

Improve this question
$\endgroup$
2
  • 1
    $\begingroup$ "numerical integrals with uncertainties" - since I am not familiar with how these are carried out in general, out of curiosity could you point to a source of the general method? $\endgroup$
    – MarcoB
    Feb 14 at 14:50
  • $\begingroup$ Since it doesn't look like Interpolation is making any use of your uncertainties and thus your integral will make no use of your uncertainties, perhaps reference.wolfram.com/language/howto/… might be of interest to you, but I still don't think your uncertainties are going to be directly applicable to NIntegrate $\endgroup$
    – Bill
    Feb 14 at 15:50

2 Answers 2

Reset to default
6
$\begingroup$

Resampling

Here I ReplaceAll Around with a RandomVariate of the NormalDistribution with the same parameters and make a list with $1000$ repeats. For each I do the Interpolation and NIntegrate obtaining a number each time. Then,

enter image description here

Around[list] gives an approximate object around the mean of the elements of list and with an uncertainty corresponding to their standard deviation.

int[y_?NumericQ] := Around@Table[
    Module[
    {
        if = Interpolation@ReplaceAll[
            data,
            Around[a_,b_] :> RandomVariate[ NormalDistribution[a,b] ]
        ]
    },
    NIntegrate[ if[x], {x, 0, y}]
    ]
    , {1000}
]

enter image description here

Now I compare with the analytical integral $y^4/4$

Show[
    ListPlot[
        Table[
            { k, int[k] }
            , {k, 0,2, 0.2}
        ]
        , PlotRange -> {0.0001, All}
        , PlotTheme -> "Scientific"
        , ScalingFunctions -> "Log"

    ],
    Plot[
        y^4/4
        , {y, 0,2}
        , ScalingFunctions -> "Log"
    ]
]

enter image description here

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you, this is precisely the solution I also came up with $\endgroup$
    – Kai
    Feb 14 at 17:02
  • 1
    $\begingroup$ One thing to add to this though, it is worthwhile to check for convergence of the mean and uncertainty as the number of resampled points is varied, although I found 1000 to be more than enough $\endgroup$
    – Kai
    Feb 14 at 17:57
5
$\begingroup$

Using the Trapezoidal_rule

Here the interpolation is happening between the two Nearest points using Rescale[x,{min,max},{ymin,ymax}].

The "Integration" is the trapezium area $(b-a) \times \left(f(a)+f(b)\right)/2$.

ClearAll[trap]; (* Trapezoidal integration *)
trap[{a_,fa_},{b_, fb_}]:=(b-a)*(fa+fb)/2

int2[y_]:= Module[
    {i,j}, (* Index of the two nearest points *)
    {i,j} = Sort@Nearest[ data[[All,1]]->"Index", y, 2];
    Rescale[
        y,
        data[[{i,j},1]],
        {#,#+trap@@data[[i;;j]]}&@Total[trap@@@Partition[data[[;;i]], 2,1]]
    ]
]

Now I compare with the analytical integral.

Show[
    ListPlot[
        Table[
            { k, int2[k] } (* The function of interest *)
            , {k, 0,2, 0.2}
        ]
        , PlotRange -> {0.0001, All}
        , PlotTheme -> "Scientific"
        , ScalingFunctions -> "Log"
        ],
    Plot[
        y^4/4  (* The analytical solution  *)
        , {y, 0,2}
        , ScalingFunctions -> "Log"
    ]
]

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.