How to compute numerical derivative after numerical integral?

Question

Here's a simple example code to demonstrate my issue:

f[x_, a_] := a x^2;
g[a_] := NIntegrate[f[x, a], {x, 0, 2}];
ND[g[a], {a, 0}]

Running this tells me that the integrand evaluates to non-numerical values. I'm assuming this is a problem of the value a being used symbolically at some point. The function I'm actually interested in integrating is much more complicated and involves integration over multiple variables - specifically, I really do want this order of computation, where Mathematica numerically evaluates the integral then takes the derivative numerically. Symbolically taking the derivative first won't work for my purposes, since the resulting integrand has some nasty singular behavior at the origin - the integral after differentiation should still technically converge, but Mathematica starts throwing a lot of warnings, so I figure it will be more numerically stable to evaluate the integral first and then take the derivative after.

Answer 1

f[x_, a_] = a x^2;
g[a_?NumericQ] := 
   NIntegrate[f[x, a], {x, 0, 2}, WorkingPrecision -> 30, 
      AccuracyGoal -> 10];

ND[g[a], a, 0, WorkingPrecision -> 33]

(*   2.6666666666666666666666666667   *)

Answer 2

This approach may or may not work with your actual integrand, but you could try using ParametricNDSolveValue instead of ND:

f[x_,a_]:=a x^2
pf = ParametricNDSolveValue[
    {g'[x] == f[x, a], g[0]==0},
    g[2],
    {x, 0, 2},
    a
];

Then:

pf'[1]

2.66667