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][[1]]
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.