# Eliminate some temporaries from a set of equations

I have a long set of equations, mostly obtained by introducing physics related quantities which stand for some kind of temporary variables:

eqns = {
lx == Log[10., en] + 0.3,
xmax == 0.849 + 0.0537 lx,
nmu == 1.37 en^-0.04,
fem == ((s - xmax + 0.689)/0.802)^(0.802/0.113) Exp[(0.113 - s + xmax)/0.113],
fmu == ((s - xmax + 0.228)/0.573)^(0.573/1.267) Exp[(0.345 - s + xmax)/1.267],
cic == 22.3 fem + 15.6 nmu fmu,
en == (s1000/cic)^(1/0.97)
};


This is 7 equations with 9 variables. I would like to eliminate all of the variables except for {en, s, s1000} and then obtain the solution in the form en[s, s1000]. I tried Solve, Reduce, Eliminate, tried taking them pairwise etc. but while there are some results for equation pairs, for the whole system Mathematica just runs for ever...

For example, while this works as expected

Eliminate[eqns[[{1, 2}]], lx]


giving 0.86511 + 0.0233216 Log[en] == xmax, trying one more

Eliminate[eqns[[{1, 2, 4}]], {lx, xmax}]


returns a result which imho does not make any sense. Looking around StackExchange one finds Mathematica (at least in v9) has a hidden 3rd argument for elimination

Solve[eqns[[{1, 2}]], xmax, lx]


which gives a warning that lx is not a proper domain specification and that lx will be used for elimination. Unfortunately, trying this trick on three equations

Solve[eqns[[{1, 2, 4}]], fem, {lx, xmax}]


returns an empty list while we would like to see a rule for fem expressed in terms of en and s. Please, help! I remember in older versions Solve used to be more flexible and more like a universal tool for solving. How particularly equation-solving cases are divided into usage of Solve, Eliminate, Reduce etc. is also not very clear, neither from the Help nor from tutorials, mostly since too simple examples are used as illustrations. A word from an expert on this subject is very welcome...

One can define a set of rules

rules = Thread[eqns[[All, 1]] -> eqns[[All, 2]]]


and then get an equation containing only en, s, s1000.

impl = eqns[] /. rules[] /. rules[] /. rules[] /.
rules[] /. rules[] /. rules[]


Now we can use FindRoot to get a solution :

sol[x_?NumericQ, y_?NumericQ] := FindRoot[impl /. {s -> x, s1000 -> y}, {en, 0.5}][[1, 2]]

Plot3D[sol[x, y], {x, 0.5, 1}, {y, 0, 1}, PlotRange -> All] • Thanks! This gives me at least an implicit expression. Nevertheless, the example above introduces all new variables in explicit expressions containing only previously defined quantities, that's why applying the rules in reverse order solves it. What if I have a system of (nonlinear?) implicit equations, is there a general way to eliminate nuisance variables? (other than trying your kind of solution manually). Tanks again! Jul 15, 2013 at 12:05
• I don't know of any generic method and I doubt it actually exists. Jul 15, 2013 at 12:09