# Incorrect numerical derivative of function that uses FindRoot

I am trying to plot the derivative of function g[x] below where g[x] is defined as the root of another equation. However, I am getting nonsensical answers for the derivative.

Any advice would be much appreciated!

Simplified code is below:

jRa[T3_] := 10^-5/(1 + 5*10^-8/T3)

jRb[x_, T3_] :=
x*1.3*10^-4/(10^-4 +
x*(1 + 10^-6/(10^-5 - T3)))

g[x_?NumericQ] := (T3 /.
FindRoot[{jRa[T3] == jRb[x, T3]}, {{T3, 10^-7}},
MaxIterations -> 1000])/Ttot

s1 = Plot[(g[x]), {x, 10^-7, 2*10^-5}, PlotRange -> All]

s2 = Plot[g'[x], {x, 10^-7, 2*10^-5}, PlotRange -> All]

g'[.9*10^-5]

• Welcome to Mathematica.SE! I have formatted the question in the way that is customary on this site. You can click the EDIT link above to see how I did this so you can do it yourself next time. Also, can you please indicate what answers you are getting yourself, and (if necessary) explain why they are wrong? It might happen that the results are not the same on others' machines, so it is good practice to do this. Also, to make the code complete, provide the numerical value of Ttot (I know that this is just a small oversight). Oct 30, 2013 at 20:41
• Now in this case, I can see the problem on my machine. The derivative Mma computes, using g'[x], is clearly incorrect and it's negative where it should be positive. It may be because of the unusual magnitude of the numbers you used, but it's clearly incorrect on my machine. The ND function from the <<NumericalCalculus  package is usually much more reliable. However, in this case it also gives a wrong result on my machine! I do not know why. Looks like a bug to me. Oct 30, 2013 at 20:43
• For other SE members, here's what I get for the function and its derivatives using various methods: dropbox.com/s/7tnzgjfjt60ykpz/derivatives.pdf (Ttot = 1.0) Oct 30, 2013 at 20:46

It is a pain to see ND behave in a erratic manner for this example! However a quick fix is also within reach! We can use finite difference here. For more detail about this check here. In the code supplied by the OP I have assumed Ttot = 1.0.

n = 1000;
h = 2*10^-5/n;
grid = h Range[0, n];
(* Interpolating function of your derivative*)
fixder[fun_, degree_, grid_] :=
Interpolation@(Transpose@{grid,
NDSolveFiniteDifferenceDerivative[Derivative[degree], grid,Map[fun, grid]]});
(* the first derivative  and the second one*)
{d1,d2} = fixder[g, #, grid]&/@{1,2};
GraphicsRow[{s1, Plot[d1[x], {x, 0, 2*10^-5}, PlotRange -> All],
Plot[d2[x], {x, 0, 2*10^-5}, PlotRange -> All]}, ImageSize -> 600] You can use implicit differentiation to find the derivative of g:

Ttot = 1.0;

dg = (Dt[T3]/Ttot)/Dt[x] /.
First@Solve[Dt[jRa[T3] == jRb[x, T3]], Dt[T3]] /. T3 -> g[x] // Simplify

(* (99.6169 (5.*10^-8 + 1. g[x])^2 (1.*10^-10 - 0.00002 g[x] +
1. g[x]^2))/(3.83142*10^-21 + 8.42912*10^-17 x +
4.66092*10^-13 x^2 + (-7.66284*10^-16 - 1.6092*10^-11 x +
1.53257*10^-8 x^2) g[x] + (3.83142*10^-11 + 7.66284*10^-7 x + 1. x^2) g[x]^2) *)

Plot[g[x], {x, 10^-7, 2*10^-5}, PlotRange -> All]
Plot[dg, {x, 10^-7, 2*10^-5}, PlotRange -> All]  