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

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.

-
Since NSolve[] will return a list of lists, you'll likely want to put Part[] in there somewhere... – J. M. Jun 17 '13 at 17:45
– Dr. belisarius Jun 17 '13 at 17:49
@bel, FWIW: it's one of the songs I listen to after a hard day... precisely that medley. – J. M. Jun 17 '13 at 17:51
Just yesterday I posted an answer where this is illustrated. There are countless other answers on this site that do similar things. Maybe just searching for NDSolve will turn up other good examples. – Jens Jun 17 '13 at 17:56
you can easily do what you asked, eg.f[a_] := x /. NSolve[{x + y == a, x^2 + y^2 == 1}, {x, y}]. I cant understand what your intent is with those examples though. The arg x on input does what? – george2079 Jun 17 '13 at 18:37

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}]

-
While this works, I would be careful about explicitely making eq a function of b: 'eq[b_] = x^2 - b*x - 1 == 0' and I would force the definition of a[b] to take only numeric values: a[b_?NumberQ]. In this example it is fine, but you have to be careful of the order of evaluation when you are mixing both symbolic and numerical values. – Jonathan Shock Jun 18 '13 at 8:41
Isn't it inefficient to call Solve every time the function a is evaluated? If you used Set instead of SetDelayed, the equation would only be solved once. – nikie Jun 18 '13 at 9:07
@nikie yes in this simple example. There are however situations where you need to Solve (or especially NSolve ) for the specific value of the parameters. – george2079 Jun 18 '13 at 19:15

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).$$

1. 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]}}

-