# How to find derivative of a numerical solution, where precision is ambiguous?

I am trying to take the derivative of a numerical solution. I am concerned that the way I'm doing this may be problematic due to numerical error; I think there must be a better way but I'm not very experienced in Mathematica. The relevant function is called GammaFcn and is defined as follows:

K = 10; c = 1; a = 1; b = 3; d = .01;  tradeoff = 1; mig = .002;
N1[x_] := K*(1 - c*x + b*x^a - mig);
f[x_, y_, t_] := 1/(1 + N1[x]*E^(c*t*(y - x)));
Gstar[x_, y_, t_] := s*(y*f[x, y, t] + x*(1 - f[x, y, t]));
sStar[x_, y_, t_] :=
s /. NSolve[(1 - s/K) + b*(Gstar[x, y, t]/s)^a -
x*c*(1 - f[x, y, t]) - y*c*f[x, y, t] - mig == 0, s][]
N2[x_, y_, t_] := f[x, y, t]*sStar[x, y, t];
Beta1[x_] := mig*N1[x];
Estar[x_] := (d)/(Beta1[x]);
Istar[x_] := 1 - Estar[x];
lambda[x_] := Beta1[x]*Istar[x];
I1star[SigmaI_, x_] :=     Estar[x]*(1 - Estar[x])*mig*N1[x]/(d + SigmaI*lambda[x]);
M1[SigmaI_, p_, q_, t_, z_] := Sum[z^(n)*(d + SigmaI*Istar[p]*mig*N1[p])* E^(-t*(d + SigmaI*Istar[p]*mig*N1[p]))* E^(-N1[q]*mig*t)*(N1[q]*mig)^n*t^n/n!, {n, 0, Infinity}]
M2[p_, q_, t_, z_] := Sum[z^(n)*((d)*E^(-t*(d))* E^(-mig*(sStar[q, p, t] - N2[q, p, t])*t)*(mig*(sStar[q, p, t] - N2[q, p, t]))^n*t^n/n!), {n, 0, Infinity}];
M3[p_, q_, t_, z_] := Sum[z^(n)*((d)*E^(-t*(d))* E^(-mig*N2[p, q, t]*t)*(mig*N2[p, q, t])^n*t^n/n!), {n, 0, Infinity}];
Part1[SigmaI_, p_, q_, z_] := NIntegrate[M1[SigmaI, p, q, t, z], {t, 0, Infinity}]
Part2[p_, q_, z_] := NIntegrate[M2[p, q, t, z], {t, 0, Infinity}]
Part3[p_, q_, z_] := NIntegrate[M3[p, q, t, z], {t, 0, Infinity}]
Const1[SigmaI_, Res_, New_] := (d/(d + SigmaI*Istar[Res]*mig*N1[Res]));
Correct[SigmaI_, Res_, New_, Z_] := Const1[SigmaI, Res, New]*Part1[SigmaI, Res, New, Z] + (1 - Const1[SigmaI, Res, New])*(Part1[SigmaI, Res, New, Z]*
Part2[Res, New, Z])

GammaFcn[SigmaI_, p_, q_, z_] := (1 - Estar[p] - I1star[SigmaI, p]) + Estar[p]*Correct[SigmaI, p, q, z] + I1star[SigmaI, p]*(SigmaI*Part3[p, q, z] + (1 - SigmaI));


I am interested in solutions to GammaFcn[SigmaI, p, q, z]=z, where 0>z>1; suppose these solutions are given by a function GammaFcnFixedPt[SigmaI, p, q]. In particular, I want to calculate the first derivative of GammaFcnFixedPt with respect to p at the point p=q, with some fixed SigmaI and q. I do this with the following function, for some small epsilon (because when p=q, GammaFcn=1):

GammaFcnFixedPtDeriv[SigmaI_, q_, epsilon_] := (-1 + z /. FindRoot[GammaFcn[SigmaI, q + epsilon, q, z] - z, {z, .99}])/  epsilon


My problem is that when I do this, when epsilon gets sufficiently small, the results diverge. So with decreasing epsilons, the function will appear to be converging onto the derivative, but then (presumably due to precision errors), it diverges. Here are pairings of epsilon and the approximated derivative (from above function), for SigmaI=.1, q=.3. :

{1/10, -0.200505}, {1/10^2, -1.44905}, {1/10^3, -0.901373}, {1/10^4, -0.886587}, {1/10^5, -0.886699}, {1/10^6, -0.886358}, {1/10^7, -0.847299}, {1/10^8, -0.297589}, {1/10^9, 0.969691}

So in the above pairings, before epsilon goes below 1/10^7, there it looks as though the sequence is approaching somewhere around .88, but thereafter diverges quickly. Again I assume this is due to some machine error; but I'm not positive it is even due to this error. But if this is the case, how do I know what epsilon to pick for the most accurate derivative?

If I add "PrecisionGoal -> 10" to the integrals in the code in the first block above:

Part2[p_, q_, z_] := NIntegrate[M2[p, q, t, z], {t, 0, Infinity}, PrecisionGoal -> 10]
Part3[p_, q_, z_] := NIntegrate[M3[p, q, t, z], {t, 0, Infinity}, PrecisionGoal -> 10]


{1/10^1, -0.200505}, {1/10^2, -1.44905}, {1/10^3, -0.901373}, {1/10^4, -0.886587}, {1/10^5, -0.886699}, {1/10^6, -0.886431}, {1/10^7, -0.855601}, {1/10^8, -0.0572474}, {1/10^9, -13.9658}

And doing the same thing with PrecisionGoal->12 gives:

{1/10^1, -0.200505}, {1/10^2, -1.44905}, {1/10^3, -0.901373}, {1/10^4, -0.886587}, {1/10^5, -0.886699}, {1/10^6, -0.886431}, {1/10^7, -0.855601}, {1/10^8, -1.83026}, {1/10^9, -13.9658}

When I set it much higher than that it takes a long time. So how can I know that the derivative is converging to (in this case) approximately -0.887, rather than convergence not occurring (limit not existing)? And is there a better way to find the derivative than I have done? Thanks very much.

You are going to get clobbered by an ugly blend of cancellation and roundoff error here. Those fixed points are on the order of epsilon away from one. If you work at machine precision then at epsilon around 10^(-8) you can expect:

8 digits at the front lost to cancellation Some number lost at the back to roundoff error. With default settings of NIntegrate this too can easily be as much as 8 digits.

So you may well get few or no correct digits in your estimated derivative.

A computationally slow remedy is as follows. Use exact inputs everywhere. Set WorkingPrecision and PrecisionGoal relatively high. And hope for the best.

In the code below I use Set rather than SetDelayed in several places, in the hope of making it a bit faster.

K = 10; c = 1; a = 1; b = 3; d = 1/100; tradeoff = 1; mig = 1/500;
N1[x_] = K*(1 - c*x + b*x^a - mig);
f[x_, y_, t_] = 1/(1 + N1[x]*E^(c*t*(y - x)));
Gstar[x_, y_, t_] = s*(y*f[x, y, t] + x*(1 - f[x, y, t]));

sStar[x_, y_, t_] =
s /. Solve[(1 - s/K) + b*(Gstar[x, y, t]/s)^a -
x*c*(1 - f[x, y, t]) - y*c*f[x, y, t] - mig == 0, s][];
N2[x_, y_, t_] = f[x, y, t]*sStar[x, y, t];
Beta1[x_] = mig*N1[x];
Estar[x_] = (d)/(Beta1[x]);
Istar[x_] = 1 - Estar[x];
lambda[x_] = Beta1[x]*Istar[x];
I1star[SigmaI_, x_] =
Estar[x]*(1 - Estar[x])*mig*N1[x]/(d + SigmaI*lambda[x]);

M1[SigmaI_, p_, q_, t_, z_] =
Sum[z^(n)*(d + SigmaI*Istar[p]*mig*N1[p])*
E^(-t*(d + SigmaI*Istar[p]*mig*N1[p]))*
E^(-N1[q]*mig*t)*(N1[q]*mig)^n*t^n/n!, {n, 0, Infinity}];
M2[p_, q_, t_, z_] :=
Sum[z^(n)*((d)*E^(-t*(d))*
E^(-mig*(sStar[q, p, t] - N2[q, p, t])*
t)*(mig*(sStar[q, p, t] - N2[q, p, t]))^n*t^n/n!), {n, 0,
Infinity}];
M3[p_, q_, t_, z_] :=
Sum[z^(n)*((d)*E^(-t*(d))*
E^(-mig*N2[p, q, t]*t)*(mig*N2[p, q, t])^n*t^n/n!), {n, 0,
Infinity}];

Part1[SigmaI_, p_, q_, z_] :=
NIntegrate[M1[SigmaI, p, q, t, z], {t, 0, Infinity},
WorkingPrecision -> 30, PrecisionGoal -> 18]
Part2[p_, q_, z_] :=
NIntegrate[M2[p, q, t, z], {t, 0, Infinity}, WorkingPrecision -> 30,
PrecisionGoal -> 18]
Part3[p_, q_, z_] :=
NIntegrate[M3[p, q, t, z], {t, 0, Infinity}, WorkingPrecision -> 30,
PrecisionGoal -> 18]

Const1[SigmaI_, Res_, New_] = (d/(d + SigmaI*Istar[Res]*mig*N1[Res]));

Correct[SigmaI_, Res_, New_, Z_] :=
Const1[SigmaI, Res, New]*
Part1[SigmaI, Res, New,
Z] + (1 - Const1[SigmaI, Res, New])*(Part1[SigmaI, Res, New, Z]*
Part2[Res, New, Z])

GammaFcn[SigmaI_?NumberQ, p_?NumberQ, q_?NumberQ,
z_?NumberQ] := (1 - Estar[p] - I1star[SigmaI, p]) +
Estar[p]*Correct[SigmaI, p, q, z] +
I1star[SigmaI, p]*(SigmaI*Part3[p, q, z] + (1 - SigmaI));

GammaFcnFixedPtDeriv[SigmaI_, q_,
epsilon_] := ((z /. (xxx =
FindRoot[GammaFcn[SigmaI, q + epsilon, q, z] - z, {z, 99/100},
WorkingPrecision -> 40, PrecisionGoal -> 15])) - 1)/epsilon


Notice I keep track of the FindRoot value; this is for later display. Here are a few representative runs.

{GammaFcnFixedPtDeriv[1/10, 3/10, 10^(-4)], N[xxx]}

(* Out= {-0.886586643003313812325995027319719947, {z ->
0.999911341336}} *)

{GammaFcnFixedPtDeriv[1/10, 3/10, 10^(-6)], N[xxx]}

(* Out= {-0.8867304424355187517516855638775763, {z -> 0.99999911327}} *)

{GammaFcnFixedPtDeriv[1/10, 3/10, 10^(-12)], xxx}

(* Out= {-0.8867337144809477198576448149, {z ->
0.9999999999991132662855190522801423551851}} *)


Among other things this may give a sense of how fast the derivative approximation converges in epsilon. It appears to be somewhat worse than linear.

The code above fares poorly for epsilon much smaller. Possibly it would improve, at the expense of speed, by further raising the option parameter values, or at least the PrecisionGoal parameter.

• Thank you very much Daniel. This seems to do the trick, and lets me confirm that the divergence experienced previously indeed are due to approximation errors. "NIntegrate::inumr: The integrand (E^(-((1103900001 t)/25000000000)-(t (<<9>>/<<7>>-<<1>>))/25000) (250 E^(304900001 t/25000000000)+(54900001 E^<<1>>)/1000000))/25000 has evaluated to non-numerical values for all sampling points in the region with boundaries {{[Infinity],0}}. >>" – mike Mar 7 '13 at 15:51
• Sorry, I meant to ask whether the warnings shown (in both my original version and your improved vs) are a problem or if I can safely ignore them: NIntegrate::inumr: "The integrand (E^(-((110390000001 t)/2500000000000)-(t (<<1>>-<<1>>))/25000)(250\E^(30490000001 t/2500000000000)+(5490000001 E^<<1>>)/100000000))/25000 has evaluated to non-numerical values for all sampling points in the region with boundaries {{[Infinity],0}}." Thanks again. – mike Mar 7 '13 at 16:22
• They may arise from symbolic preprocessing, in which case they are benign. I don't think I was getting them though because I had defined a restricted GammaFcn that only fired for explicit NumberQ arguments. – Daniel Lichtblau Mar 7 '13 at 17:44