How could I define a function as the solution of equations in Mathematica?

Question

I have posed a specific question yesterday but it may be too trivial to answer.

I think the question could be asked in a more general way so it is easier to answer. Then I could solve the original problem by myself.

Usually, we define a function as

F[x_] := x + 10 

where the right-hand side is an expression.

Now, I would like define a function where the right-hand side is the solution of a set of equations. Like

optV[V_, r_, λ_, μ_, η_, σ_] := 
  V /. 
    FindRoot[{
      A V^θ[r, λ, μ, η, σ] h[V, c, r, λ, μ, η, σ] == V - c, 
      D[A*V^θ*h[V, c, r, λ, μ, η, σ], V] == 1
      }, 
      {A, 0}, {V, 1}]; 

a[V_, r_, λ_, μ_, η_, σ_] := 
  A /. 
    FindRoot[{
      A V^θ[r, λ, μ, η, σ] h[c, r, λ, μ, η, σ] == V - c, 
      D[A V^θ[r, λ, μ, η, σ]  h[V, c, r, λ, μ, η, σ], V] == 1
      }, 
      {A, 0}, {V,1}]

This may not be correct code, but I think you can understand my question from it. optV and a are both parameters for another final function which would be manipulated. So you could regard optV and a as intermediate parameters. So the solution of the set of equations, correspondingly optV and a, varies along with c, r, λ, μ, η, σ, which are parameters input for manipulating plot.

Answer 1

If I understand you right, here is the answer, but it is in a form of a more transparent example.

Let this

Clear[a,b];
eq=x^2 - b*x - 1 == 0

be an equation depending upon a parameter b with b>0. Let us define a=a[b] which is its smaller solution of the equation eq:

a[b_] := Solve[eq, x][[1, 1, 2]]

It seems that this (or something alike) is what you need. You can check that a is indeed a function of b by, say, plotting a[b]. Evaluate this:

Plot[a[b], {b, 0, 3}]

With some care one can do the same with the FindRoot statement:

Clear[a,b];
a[b_] := FindRoot[eq, {x, 0}][[1, 2]]
Plot[a[b], {b, 0, 3}]

Answer 2

I will give an elaborate answer because I think this is an area where Mathematica really shines ;)

Let me digress.

We know that, in general, an equation system defines a set of points as a function of the parameters in it. Namely, the set of variable assignments or valuations that make the equation true.

Thus an equation system can also be seen as a predicate $E(\vec x, \vec p)$ which defines the set of its solutions given the parameters $\vec p$: $$sol(E, \vec p) := \{\vec x : E(\vec x, \vec p)\}$$

Note that "solving" an equation system, even for given parameter values (for example if you have no parameters) does not usually give a single value $\vec x$.

Now what does mathematica's Solve[E, {x1,...,xn}] (or NSolve) do? It gives a representation of this set in the form of a finite list of replacement rule sets, $\{r_1, r_2, ..., r_n\}$.

Each $r_i$ contains rules of the form $x_{ij} \mapsto f_{ij}(\vec x', \vec p)$. Let $(\vec x)_i'$ be the sublist of $\vec x$ not appearing as an $x_{ij}$ (on the left hand side of a rule), $(\vec x)_i''$ be the others. Then the ruleset $r_i$ defines a map $F_i$ with $$(\vec x)_i'' := F_i((\vec x)_i', \vec p).$$ The set of solutions is then characterized as follows: $$\vec x \in sol(E, \vec p) \iff \exists i\quad (\vec x)_i'' = F_i((\vec x)_i', \vec p).$$

---

Turning this into a function

Ideally, we would want to get a function $G$ with

$$G(\vec p) := sol(E, \vec p)$$

Obviously, this is in general impossible because the set of solutions is infinite (sometimes not even countable).

But assume we knew $E$ has a finite amount of solutions for every fixed $\vec p$. Can we derive these from the output of Solve?

In general no, one because of the parameters $\vec C$ mathematica might add to describe the solution space and because we dont know which $(\vec x)_i'$ will not generate Undefined values with the functions we are given. And we get no explicit way of enumerating only valid candidates.

Still, if the $F_i$ are total in $(\vec x)_i'$ we might be able to leverage them to describe our solution space in a more explicit way.

In particular, we can ask Mathematica to give us the solution in 1 variable in terms of the parameters and the others, simply by 'solving' for that variable only. Using ParametricPlot and friends you can then visualize the solution.

---

Wrap up

Lets assume you know that $(\vec x)_i'' = \vec x$ for all $r_i$, i.e. that all variables you care about have been explicitly solved for. Then the function $$f(\vec p) := \{F_1(\vec p), ...,F_n(\vec p)\}$$ can be defined in the following two ways in mathematica. (Let us write {{x1, ...}*} = f[{p1, p2,...}]).

1. Solving the equation system once

If you know (or enforce) that x1,...,p1,... are undefined at the point of definition:

ClearAll[f, p1, p2, x1, x2];
(*Define f[{p1,...,pn}]. {x1,...,xn, p1,...,pn} must not be defined at this point \
for this to work*)
f[{p1_, p2_}] := Evaluate[{x1, x2} /. Solve[
     (*The equation system*)
     x2 == p1^2 + p1 && p2 == x1,
     {x1, x2}]];
(*Test*)
?f
f[{3, 4}]

Otherwise (namespace-clean version):

ClearAll[f];
(*Prepare f*)
Evaluate[Module[{x1, x2, p1, p2},
   f[params] = {p1, p2};
   f[variables] = {x1, x2};
   f[solutions] = f[variables] /.
     Solve[
      (*The equation system*)
      x2 == p1^2 + p1 && p2 == x1, f[variables]]]
  ];
(*Define f[{p1,...,pn}]*)
f[p_] := f[solutions] /. Rule @@@ Transpose@{f[params], p};
(*Test*)
?f
f[{3, 4}]

2. Solving the equation system on every call

Namespace-dirty version -- x1,... must not be defined when calling f.

ClearAll[f, x1, x2];
(*When calling this function, x1,...,xn must be undefined*)
f[{p1_, p2_}] := {x1, x2} /. Solve[
    (*The equation system*)
    x2 == p1^2 + p1 && p2 == x1,
    {x1, x2}];
(*Test*)
?f
f[{3, 4}]

To make this clean, simply wrap the definition in a module:

ClearAll[f, x1, x2];
Module[{x1, x2},
  f[{p1_, p2_}] := {x1, x2} /. Solve[
      (*The equation system*)
      x2 == p1^2 + p1 && p2 == x1,
      {x1, x2}];
  ];
(*Test*)
?f
f[{3, 4}]

---

Some more details

1. Note that the $F_i$ are in general partial functions. You may get

   x -> ConditionalExpression[1, a > 0]

   which is Undefined for certain values, try it:

   ConditionalExpression[1, a > 0] /. {a -> 0}

2. The above explanation makes it clear why {} denotes the empty set of solutions.

   And also why {{}} is the full dimensional set of solutions: In this case $(\vec x)_1' = \vec x$, $(\vec x)_1''$ is empty and the function $F_1$ returns nothing, thus $(\vec x)_1'' = F_1((\vec x)_1',\vec p)$ for any $\vec x$.

3. Sometimes Solve will introduce new parameters $\vec C$ in the solution. These parametrize the $F_i$ so that actually, implicitly, infinitely many $F_i$ are available to pick to construct a solution:

$$\vec x \in sol(E, \vec p) \iff \exists i\exists \vec C\quad (\vec x)_i'' = F_i((\vec x)_i', \vec p, \vec C).$$

4. Because "Solve gives generic solutions only. Solutions that are valid only when continuous parameters satisfy equations are removed" the picture I painted above is a bit too simple. Sometimes there will be solutions that Mathematica is well aware of but does not give because of this. It will generate ConditionalExpression with inequalities for parameters:

   Solve[x^2 == a - b, x, Reals]

   {{x -> ConditionalExpression[-Sqrt[a - b], a > b]}, {x -> ConditionalExpression[Sqrt[a - b], a > b]}}

   but by default 'forgets' all solutions that would require an equality between parameters:

   Solve[x == 0 && x^2 == a - b, x, Reals]

   {}

   It will also generate solutions that are wrong for a finite set of parameter assignments: Solve[x a == 1 , {x, y}] {{x -> 1/a}}

   I don't see why you shouldn't always use MaxExtraConditions -> All which gives these solutions back:

   Solve[x == 0 && x^2 == a - b, x, Reals, MaxExtraConditions -> All]

   {{x -> ConditionalExpression[0, a == b]}}

   Solve[x a == 1 , {x, y}, MaxExtraConditions -> All]

   {{x -> ConditionalExpression[1/a, a != 0]}}