Total derivative after numerical solution of a system

I have a system of three equations in three unknowns, $$k$$, $$\theta$$ and $$w$$. I am interested in the behavior of the variables when one parameter changes. I first specify the system:

Clear[eq1, eq2, eq3, k, theta, w, a, eta, r, s, p, rho, beta, solk,
soltheta, solw];

eta = 0.5; a = 1/3; r = 1/20; s = 1/10; p = 1/5; rho = 0.7;

f = (a*k^rho + (1 - a))^(1/rho);

eq1[beta_] :=
p - ((1 - beta)*theta^(-eta)/(r + s + (1 - beta)*theta^(-eta))*
D[f, k]/(r + s));

eq2[beta_] :=
p - ((1 - beta)*
theta^(-eta)/(r + s +
beta*theta^(1 - eta) + (1 - beta)*theta^(-eta))*f/k*1/(r + s));

eq3[beta_] :=
w - (beta*(
r + s + theta^(1 - eta))*(a*k^rho + (1 - a))^(1/rho)/(r + s +
beta*theta^(1 - eta) + (1 - beta)*theta^(-eta)));

Then, I solve the system numerically and I plot the solution as functions on the parameter $$\beta$$:

solk[beta_?NumericQ] := {k} /.
FindRoot[{eq1[beta] == 0, eq2[beta] == 0,
eq3[beta] == 0}, {{k, 0.00009}, {theta, 1000000}, {w, 0.6}}];

soltheta[beta_?NumericQ] := {theta} /.
FindRoot[{eq1[beta] == 0, eq2[beta] == 0,
eq3[beta] == 0}, {{k, 0.00009}, {theta, 1000000}, {w, 0.6}}];

solw[beta_?NumericQ] := {w} /.
FindRoot[{eq1[beta] == 0, eq2[beta] == 0,
eq3[beta] == 0}, {{k, 0.00009}, {theta, 1000000}, {w, 0.6}}];

Plot[solk[beta], {beta, 0.1, 0.9}]
Plot[soltheta[beta], {beta, 0.1, 0.9}]
Plot[solw[beta], {beta, 0.1, 0.9}]

Now, I would like to get some insight about why, in particular, $$w$$ behaves like it does.

For instance, I could compute the total derivative of eq3:

Dt[eq3[beta]] // FullSimplify

But the output is rather complicated and difficult to interpret. What I would like to obtain is an expression in which $$\frac{\partial{w}}{\partial{\beta}}$$ depends on $$\frac{\partial{k}}{\partial{\beta}}$$ and $$\frac{\partial{\theta}}{\partial{\beta}}$$. How can I do that?

• Do you need dependencies of derivatives on solutions of equations 1,2,3? – Alex Trounev Nov 5 '18 at 10:41
• Maybe D[w /. First@Solve[eq3[beta] == 0, w] /. {k -> k[beta], theta -> theta[beta]}, beta] // Simplify? You could use Rationalize[eq3[beta]] to simplify the numbers, which seem to be simple fractions. Or if you like Dt, then perhaps Dt[w /. First@Solve[Rationalize@eq3[beta] == 0, w]] /. {Dt[beta] -> 1} // Simplify? – Michael E2 Nov 5 '18 at 12:48
• @Alex I guess so. Since the system is block recursive I would like to know how the derivative of $w$ wrt $\beta$ depends on $k^*(\beta)$ and $\theta^*(\beta)$ where stars means the functions solves the system eq1 and eq2 – Giop Nov 6 '18 at 12:57
• @MichaelE2, thank you. So the output of D[w /. First@Solve[eq3[beta] == 0, w] /. {k -> k[beta], theta -> theta[beta]}, beta] // Simplify gives me $w$ as a function of $k(\beta)$ and $\theta(\beta)$ and their derivatives. Unfortunately the expression is still too complicated to be studied analytically – Giop Nov 6 '18 at 14:15

It is necessary to replace the variables in the equations for solutions, and then calculate the derivative, for example,

e3 =
eq3[beta] /. {w -> solw[beta], k -> solk[beta],
theta -> soltheta[beta]}

Out[]= -((
beta (2/3 + solk[beta]^0.7/3)^1.42857 (3/20 + soltheta[beta]^0.5))/(
3/20 + (1 - beta)/soltheta[beta]^0.5 + beta soltheta[beta]^0.5)) +
solw[beta]
e3D = D[e3, beta]

Out[]= -(((2/3 + solk[beta]^0.7/3)^1.42857 (3/20 +
soltheta[beta]^0.5))/(
3/20 + (1 - beta)/soltheta[beta]^0.5 + beta soltheta[beta]^0.5)) - (
0.333333 beta (2/3 + solk[beta]^0.7/3)^0.428571 (3/20 +
soltheta[beta]^0.5) Derivative[solk][beta])/(
solk[beta]^0.3 (3/20 + (1 - beta)/soltheta[beta]^0.5 +
beta soltheta[beta]^0.5)) - (
0.5 beta (2/3 + solk[beta]^0.7/3)^1.42857 Derivative[soltheta][
beta])/((3/20 + (1 - beta)/soltheta[beta]^0.5 +
beta soltheta[beta]^0.5) soltheta[beta]^0.5) + (
beta (2/3 + solk[beta]^0.7/3)^1.42857 (3/20 +
soltheta[beta]^0.5) (-(1/soltheta[beta]^0.5) +
soltheta[beta]^0.5 - (
0.5 (1 - beta) Derivative[soltheta][beta])/
soltheta[beta]^1.5 + (0.5 beta Derivative[soltheta][beta])/
soltheta[beta]^0.5))/(3/20 + (1 - beta)/soltheta[beta]^0.5 +
beta soltheta[beta]^0.5)^2 + Derivative[solw][beta]

Other equations are calculated similarly.

e1 = eq1[beta] /. {w -> solw[beta], k -> solk[beta],
theta -> soltheta[beta]}; e1D = D[e1, beta]
e2 =
eq2[beta] /. {w -> solw[beta], k -> solk[beta],
theta -> soltheta[beta]}; e2D = D[e2, beta]