# Can't NSolve system of exponential and linear equations in several variables

At the moment I am stuck with the following code, when specified to use root , quadratic or log utility functions all it fine, only exponentials like:

$\small \left\{-e^{-c_{0,1}}-0.25 e^{-c_{1,1}}-0.25 e^{-c_{2,1}}-0.5 e^{-c_{3,1}}\right\}$ s.t $\small\left\{\text{pr}_1 \left(c_{1,1}-5\right)+\text{pr}_2 \left(c_{2,1}-10\right)+\text{pr}_3 \left(c_{3,1}-10\right)+c_{0,1}-10\right\}$

won't give out any solution, even when all equations (market-clearing like sum of agents consumption equal to endowments of state) are specified.

The whole file can be downloaded here. Every section is explained. The utility-function is modified from-to:

    Subscript[Vf, 1][j_, i_, \[Alpha]_] := consall[[i, j]] - \[Alpha]*consall[[i, j]]^2
Subscript[Vf, 1][j_, i_, \[Alpha]_] := -Exp[-consall[[i, j]]]


First one gives output of the two solutions:

$$\small \left\{\left\{c_{-1,1}\to 0.816667,c_{-1,2}\to 0.883333,c_{0,1}\to 9.16667,c_{0,2}\to 5.83333,c_{1,1}\to 9.16667,c_{1,2}\to 5.83333,c_{2,1}\to 9.16667,c_{2,2}\to 5.83333,c_{3,1}\to 50.,c_{3,2}\to 50.,\text{pr}_1\to 0.25,\text{pr}_2\to 0.25,\text{pr}_3\to 0.\right\},\left\{c_{-1,1}\to 0.824265,c_{-1,2}\to 0.875735,c_{0,1}\to 8.78673,c_{0,2}\to 6.21327,c_{1,1}\to 8.78673,c_{1,2}\to 6.21327,c_{2,1}\to 8.78673,c_{2,2}\to 6.21327,c_{3,1}\to 11.211,c_{3,2}\to 8.78896,\text{pr}_1\to 0.25,\text{pr}_2\to 0.25,\text{pr}_3\to 0.470588\right\}\right\}$$

Exponential only yield the message: "Solve was unable to solve the system with inexact coefficients or the system obtained by direct rationalization of inexact numbers present in the system. Since many of the methods used by Solve require exact input, providing Solve with an exact version of the system may help."

All parameters are given numerical(random or specified) values so that it should solve the system. FindRoot is not suitable as the number of solutions are too big to find the right one even if ranges are specified (ranges of mx(n-1)+n-1 variables would have to be stated !!) Any help would be great, Thanks!

• NSolve[] is really only meant for algebraic equations; once you throw in transcendentals, you use FindRoot[] along with a good guess, or give up in despair. – J. M. is in limbo May 4 '15 at 14:56
• FindRoot is not suitable as the number of solutions are too big to find the right one even if ranges are specified (ranges of m*(n-1)+n-1 variables would have to be stated !!) Solving by hand by the way gives unique solutions. Any further help would be great, Thanks! – chrisoutwright May 4 '15 at 14:58
• Then NSolve[] is even less suitable here. Do you really not have any good guesses for the parameters you're solving? – J. M. is in limbo May 4 '15 at 15:01
• Additionally, I see that some of your unknowns are within a decaying exponential; these things are notoriously badly behaved with respect to solution methods. – J. M. is in limbo May 4 '15 at 15:02
• I do know the solutions already for the problem in n=3 (first s present period) therefore 2 states. But stating ranges in great proximity to the solutions is not acutally what I want as most of the time I can't refer to solutions ;-( Modify this: 'b1 = {{0.5, 0.5}, {0.75, 0.25}}' and 'kalla = {{10, 5, 10}, {10, 10, 10}}' for n=3 – chrisoutwright May 4 '15 at 15:04

I'm not going to give a full answer because it will appear as a mess of Subscript notation (if you want to avoid such things, use a different notation at least until you have a debugged code base). There are at least a couple of routes one might take to get solutions. One, which I did with success, is to take the solutions to the original system as inputs to FindRoot for the revised system. From this I got two distinct solutions.

tester = vartab /. Solve[Flatten[fefe2], vartab];

eqns3

(* Out[346]= {E^-Subscript[c, 0, 1] == Subscript[c, -1, 1],
0.25 E^-Subscript[c, 1, 1] == Subscript[pr, 1] Subscript[c, -1, 1],
0.25 E^-Subscript[c, 2, 1] == Subscript[pr, 2] Subscript[c, -1, 1],
0.5 E^-Subscript[c, 3, 1] ==
Subscript[pr, 3] Subscript[c, -1, 1], -10 + Subscript[c, 0, 1] +
Subscript[pr, 1] (-5 + Subscript[c, 1, 1]) +
Subscript[pr, 2] (-10 + Subscript[c, 2, 1]) +
Subscript[pr, 3] (-10 + Subscript[c, 3, 1]) == 0,
E^-Subscript[c, 0, 2] == Subscript[c, -1, 2],
0.25 E^-Subscript[c, 1, 2] == Subscript[pr, 1] Subscript[c, -1, 2],
0.25 E^-Subscript[c, 2, 2] == Subscript[pr, 2] Subscript[c, -1, 2],
0.5 E^-Subscript[c, 3, 2] ==
Subscript[pr, 3] Subscript[c, -1, 2], -5 + Subscript[c, 0, 2] +
Subscript[pr, 1] (-10 + Subscript[c, 1, 2]) +
Subscript[pr, 2] (-5 + Subscript[c, 2, 2]) +
Subscript[pr, 3] (-10 + Subscript[c, 3, 2]) ==
0, -15 + Subscript[c, 0, 1] + Subscript[c, 0, 2] ==
0, -15 + Subscript[c, 1, 1] + Subscript[c, 1, 2] ==
0, -15 + Subscript[c, 2, 1] + Subscript[c, 2, 2] == 0} *)

start1 = Transpose[{vartab, tester[[1]]}]
rt1 = FindRoot[eqns3, start1, MaxIterations -> 1000]

(* Out[342]= {{Subscript[c, -1, 1], 0.816666666667}, {Subscript[c, -1,
2], 0.883333333333}, {Subscript[c, 0, 1],
9.16666666667}, {Subscript[c, 0, 2], 5.83333333333}, {Subscript[c,
1, 1], 9.16666666667}, {Subscript[c, 1, 2],
5.83333333333}, {Subscript[c, 2, 1], 9.16666666667}, {Subscript[c,
2, 2], 5.83333333333}, {Subscript[c, 3, 1], 50.}, {Subscript[c, 3,
2], 50.}, {Subscript[pr, 1], 0.25}, {Subscript[pr, 2],
0.25}, {Subscript[pr, 3], 0.}} *)

(* Out[343]= {Subscript[c, -1, 1] -> 0.000104464143832,
Subscript[c, -1, 2] -> 0.00292829969482,
Subscript[c, 0, 1] -> 9.16666666667,
Subscript[c, 0, 2] -> 5.83333333333,
Subscript[c, 1, 1] -> 9.16666666667,
Subscript[c, 1, 2] -> 5.83333333333,
Subscript[c, 2, 1] -> 9.16666666667,
Subscript[c, 2, 2] -> 5.83333333333, Subscript[c, 3, 1] -> 50.,
Subscript[c, 3, 2] -> 50., Subscript[pr, 1] -> 0.25,
Subscript[pr, 2] -> 0.25, Subscript[pr, 3] -> 6.56496072263*10^-54} *)


Similar for the other one.

The second possibility is to separate into polynomial and nonpolynomial equations, solve for an appropriate set of variables in the polynomial subsystem, and substitute those into the remaining equations. One way to go about this is indicated below.

eqns4 = Select[eqns3, FreeQ[#, E] &];
solns = Solve[
eqns4, {Subscript[c, 0, 1], Subscript[pr, 2], Subscript[c, 2, 2],
Subscript[c, 1, 2], Subscript[c, 3, 2]}]
eqns5 = Together[eqns3 /. solns] /. True -> Sequence[]
Solve[eqns5[[1]]]


It seems to be quite slow at the last step though and I do not know if it will complete in reasonable time.

• Hello Daniel, thanks very much for the answer! I've noticed that I used Exp[] instead of E[], which caused problems right at the beginning. Thanks for indicating. How is " eqns3" defined? – chrisoutwright May 4 '15 at 17:00
• I created eqn3 somewhat laboriously by replicating your code with the second variant of the Vf definition, that is, Subscript[Vf, 1][j_, i_, \[Alpha]_] := -Exp[-consall[[i, j]]] – Daniel Lichtblau May 4 '15 at 17:12
• I see, might it possible to have a link to the revised file? I'm a bit at a loss at how to implement your part by myself. Your idea is great though. The purpose of the code is to verify my solutions for exam practice, solutions to exp functions is therefore a great help. It is true that the subscript notation is a bit misguided way to use mathematica but the algebra looks more intuitive afterwards when output is presented. – chrisoutwright May 4 '15 at 17:21
• I'm sorry, I no longer have it. What system do you get when you use that definition for Subscript[Vf, 1]? Chances are it will be the same as my eqns3, – Daniel Lichtblau May 4 '15 at 17:23
• Actually Subscript[Vf, 1] is actually a dummy-element for generating the whole utility function which is composed of the dummy times probability. It can't be used singly as input. But I will try it. – chrisoutwright May 4 '15 at 17:31