# Finding numerical/analytic solution to a derivative

Consider the following system of equations, with 11 dependent variables:

seq={
f - h (1*^8 r + 200 + 10) - (25. g h)/(g + 244 h) == 0,
(25. g h)/(g + 244 h) - (3.7*^-5 gs)/(4.9*^-5 + gs) == 0,
g == 0.0032 - 2 gs,
rt == 701 r,
oht == 701 oh,
sst == 701 ss,
h r == 2.1*^-3 sst tr,
oht == 3.85*^8 h r,
701 (oh + r + ss) == 5.7*^-4,
sst tr == (0.41*^-6 to)/(1.48*^-6 + to),
to == 6*^-7 - tr
};


I am trying to find the numerical value of the derivative (dh/df)*(f/h) around a certain value of f = 5*^-9.

To do so, I could first find the numerical solution to the equations using FindRoot (fairly easy to do), and then replacing the numerical value of variables g and r in equation

f - h (1*^8 r + 200 + 10) - (25. g h)/(g + 244 h) == 0


while leaving h and f unevaluated and then compute the derivative. Yet, because those other variables (r and g) also change with f, I can't compute this just numerically like this.

Yet, finding the analytic solution to the equations above is not feasible as the system is too big for mathematica to find such solution (the sample above is just a small part of the real system). So, I was thinking, whether there would be a way for me to compute the derivative at the top numerically.

I am probably just complicating things, and most likely there is a rather trivial way of solving this - yet, I am not seeing it.

Could you guys help me out please?

-
Sosi -- can you make up a 2 or 3 parameter version of this that retains the essence of the problem but strips away irrelevant details? –  bill s Mar 28 '13 at 7:38
Unfortunately no. Were you thinking about normalizing and rescaling everything? –  Sosi Mar 28 '13 at 9:48
I was thinking about trying to put your Mathematica problem into a comprehensible form where someone here can help you. –  bill s Mar 28 '13 at 10:14
ok, I edited everything a bit, maybe it clarifies things. Thanks @bill s –  Sosi Mar 28 '13 at 10:28

We can look at seq as a system $F({\bf x}, t) = 0$ of a function $F \colon {\bf R}^{11+1} \rightarrow {\bf R}^{11}$ that defines ${\bf x}$ as a function of a parameter $t$. Here $F$ is the vector of differences of the two sides of each equation in seq, and $({\bf x}$, t) consists of the variables in seq:

vars = Union@Cases[seq, _Symbol, Infinity]
(* {f, g, gs, h, oh, oht, r, rt, ss, sst, to, tr} *)


Specifically for this problem, the parameter $t$ corresponds to the variable f and ${\bf x}$ consists of all the other variables, such as h. The problem is to find the derivative of h as a function of f (times f/h).

To find the derivatives of the components of ${\bf x}$, one applies the Implicit Function Theorem: Differentiate the system of equations treating the components of ${\bf x}$ as functions of $t$ and solve for the derivatives of the components of ${\bf x}$.

Below, the function $F$ is represented by fnX = Subtract @@@ seq and ${\bf x}$ by varsX. The derivative of $F$ with respect to ${\bf x}$ is represent by the matrix dseq.

seq = {f - h (1*^8 r + 200 + 10) - (25. g h)/(g + 244 h) == 0,
(25. g h)/(g + 244 h) - (3.7*^-5 gs)/(4.9*^-5 + gs) == 0,
g == 0.0032 - 2 gs, rt == 701 r, oht == 701 oh, sst == 701 ss,
h r == 2.1*^-3 sst tr, oht == 3.85*^8 h r,
701 (oh + r + ss) == 5.7*^-4,
sst tr == (0.41*^-6 to)/(1.48*^-6 + to), to == 6*^-7 - tr};
fnX = Subtract @@@ seq;  (* the function F *)
varsX = Cases[Union@Cases[seq, _Symbol, Infinity], Except[f]]; (* its variables x *)
dseq = D[fnX, {varsX}];  (* the matrix derivative of F *)


You indicated an interest in f == 5.*^-9. Mathematica finds eight solutions.

f0 = 5.*^-9; (* a particular value *)
solX0 = NSolve[seq /. f -> f0, varsX];
(* eight solutions omitted *)


Next set up the derivative of the system with respect to f and solve it for the derivatives of the variables in varsX. By the Chain Rule the derivative of fnX equals dseq.(#' & /@ varsX) + D[fnX, f]. For instance h' represents dh/df.

solDF = Solve[dseq.(#' & /@ varsX) + D[fnX, f] == 0, #' & /@ varsX];
(h' /. First@solDF) (f/h) /. f -> f0 /. solX0 // Sort
(* {-0.000129206, -0.0000101884, -6.84803*10^-6, -1.72425*10^-7,
-1.04181*10^-7, -1.03437*10^-7, -7.59983*10^-10, Indeterminate} *)


One can also use Cramer's Rule, which in this case is not too inefficient.

jac = Det[dseq]; (* Jacobian *)
hPos = Position[varsX, h][[1, 1]];
dhdf = Module[{mat = dseq},
mat[[All, hPos]] = -D[fnX, f];
Det@mat];

(dhdf/jac) (f/h) /. f -> f0 /. NSolve[seq /. f -> f0, varsX] // Sort
(* {-0.000129206, -0.0000101884, -6.84803*10^-6, -1.72425*10^-7,
-1.04181*10^-7, -1.03437*10^-7, -7.59983*10^-10, ComplexInfinity} *)

-
don't you mean "By the Chain Rule the derivative of fnX equals dseq.(#' & /@ varsX)"? –  Sosi Jul 3 '13 at 10:01
also, in solDF don't you mean solDF = Solve[dseq.(#' & /@ varsX) - D[fnX, f] == 0, #' & /@ varsX]; (attention to the minus) –  Sosi Jul 3 '13 at 10:12
@Sosi Given $F[x, y, z, ..., f]$ and $x = x[f]$, $y = y[f]$, etc, then $${dF\over df} = {\partial F\over\partial x} {dx\over df} + {\partial F\over\partial y} {dy\over df} + \cdots + {\partial F\over\partial f}\,.$$ Here $F$ is fnX, and in the previous line I meant to make it clear that the derivative is the ordinary deriv. of fnX with resp. to f. That's also why it is + D[fnX, f]. Does that make sense? –  Michael E2 Jul 3 '13 at 13:29
Yes, it does! Thank you so much! –  Sosi Jul 3 '13 at 14:17