# Applying Numerical Differentiation on the solution of a FindRoot problem

I would like to numerically differentiate the solution of a FindRoot problem in Mathematica, I have been struggling with this for two days and cannot find the way to do it. I simplify the problem below (so indeed here FindRoot is not really needed):

Suppose I have two equations:

2 b x - 3 b + 2 y = 0

b y - x - a = 0


where $$x$$ and $$y$$ are the variables to be found, and $$a$$ and $$b$$ are exogenous parameters.

Then, x = (3b^2-2a) / (2*b^2+2) and dx/da= -1 / (b^2+1)

Thus, dx/da evaluated at a=2 & b=2 is equal to -0.2.

In Mathematica the code would look something like this but it does not work to me:

Needs["NumericalCalculus"]

solution[a_?NumericQ, b_?NumericQ] :=
FindRoot[{2bx - 3b + 2y == 0, b*y - x - a == 0}, {{x, 2.}, {y, 2.}}]

ND[solution[2, 2][[1]], {x, 1}, 1]


Here is one solution:

solution[a_?NumericQ,b_?NumericQ]:=FindRoot[{2b*x-3b+2y==0,b*y-x-a==0},
{{x,2.},{y,2.}}];
xsolution[a_?NumericQ,b_?NumericQ]:=x/.solution[a,b];
ysolution[a_?NumericQ,b_?NumericQ]:=y/.solution[a,b];

(* numerical differentiation *)
NumericalCalculusND[xsolution[a,2],
{a,1}, (* <--- first derivative *)
2]     (* <--- at the point 2 *)
(* -0.2 *)


• In OPs code, ND[solution[2,2][[1]],...] note that solution[2,2][[1]] is a constant and the derivative will be zero.
• Just writing ND[solution[a,2][[1]],...] is problematic since solution[a,2][[1]] does not evaluate yet and will give the First of the expression which is a. Therefore I used xsolution to delay that evaluation.
• Also, OPs code has ND[...,{x,1},...] when what we actually want is differentiation with respect to a meaning ND[...,{a,1},...].
• great! thank you very much. This is exactly what I was looking for! very helpful! Commented Aug 2, 2022 at 18:58

One idea is to just augment your FindRoot call with information about the derivative. For your example, this might look like:

f[x_, y_] := 2 b x - 3 b + 2 y
g[x_, y_] := b y - x - a

FindRoot[
{
f[x[a], y[a]] == 0,
g[x[a], y[a]] == 0,
D[f[x[a], y[a]], a] == 0,
D[g[x[a], y[a]], a] == 0
} /. {a->2, b->2},
{x[2], 1}, {y[2], 1}, {x'[2], 1}, {y'[2], 1}
]


{x[2] -> 0.8, y[2] -> 1.4, x'[2] -> -0.2, y'[2] -> 0.4}

• thank you for the suggestion, but I needed to numerically differentiate the solutions. Commented Aug 2, 2022 at 19:00
• @my2cents It is quite often possible to avoid the need for numerical differentiation, which is good because obtaining a good result numerically can be quite finicky. Commented Aug 2, 2022 at 20:17
• Your are certainly right, but unfortunately this time I cannot avoid numerical differentiation. Commented Aug 2, 2022 at 21:51