# Assign FindRoot output vector to a vector instead of element by element

Using FindRoot to solve a multivariate equation system, assigning the output to variables one by one works, but assigning the output vector to a vector gives the error (MWE below) SetDelayed::shape: Lists {xsola,ysola} and {x,y}/. FindRoot[{D[(x-1)^2]==0,D[(y-2)^2]==0},{{x,0.3},{y,0.4}}] are not the same shape. The output of FindRoot is {1., 2.} and the list {xsola,ysola} prints as {xsola,ysola}. To me, these lists do look the same shape, but seemingly not to Mathematica.

How to assign a vector output of FindRoot to a vector?

MWE:

Clear[x, y, xsol, ysol, xsola, ysola]
xsol := x /.
FindRoot[{D[(x - 1)^2] == 0,
D[(y - 2)^2] == 0}, {{x, 0.3}, {y, 0.4}}]
ysol := y /.
FindRoot[{D[(x - 1)^2] == 0,
D[(y - 2)^2] == 0}, {{x, 0.3}, {y, 0.4}}]
{xsola, ysola} := {x, y} /.
FindRoot[{D[(x - 1)^2] == 0,
D[(y - 2)^2] == 0}, {{x, 0.3}, {y, 0.4}}]
{x, y} /.
FindRoot[{D[(x - 1)^2] == 0, D[(y - 2)^2] == 0}, {{x, 0.3}, {y, 0.4}}]
{xsola, ysola}
{xsol, ysol, xsola, ysola}


The reason to assign the output of FindRoot is to use it in fixed point calculations, similarly to this question: wrong use of _?NumberQ? The equations represented in the MWE as D[(x-1)^2]==0 etc will be solutions of integrals, often not in closed form, depending on parameters. Integrate has trouble with them and NIntegrate does not work because of the non-numeric parameters.

A related question More issues Integrate, NIntegrate, FindRoot suggests using == under FindRoot, which I am already doing.

More background in case anyone is curious and has too much time on their hands: a bigger piece of code (possibly with errors) into which I am trying to fit the FindRoot is

Clear[d, s, vi, vj, vlb, mui, muj, pi, pj, profit, ci, cj, pbr, \
pisol, pjsol, pis, pjs]
$Assumptions = 0 < {s, pi, pj, pis, pjs} < 1/2 && 0 < {mui, muj, vlb, vi, vj} < 1 && 0 <= {ci, cj} < 1; vlb[pi_, pj_] = 1 + pi - pj - Sqrt Sqrt[s]; ci = 0.; cj = 0.05; mui = 0.5; muj = 0.5; d[pi_, pj_, s_, pis_] = mui*(1 - Max[pi, vlb[pi, pj]]) + mui*Integrate[vi - pi + pj, {vi, Max[pi, 0], Max[pi, vlb[pi, pj]]}, Assumptions -> 0 < {s, pi, pj, pis, pjs} < 1/2 && 0 < {vlb} < 1] + muj*Integrate[(1 - Max[pi, vj - pj + pi]), {vj, 0, vlb[pj, pis]}, Assumptions -> 0 < {s, pi, pj, pis, pjs} < 1/2 && 0 < {vlb} < 1]; profit[pi_, pj_, s_, pis_, c_] = (pi - c)*d[pi, pj, s, pis]; pbr[pj_Real?Positive, s_Real?Positive, pis_Real?Positive, c_Real?Positive] = NArgMax[{profit[pi, pj, s, pis, c], c < pi < 1/2}, pi] pisol[s_Real?Positive] := pi /. FindRoot[{pbr[pj, s, pis, ci] == pi, pbr[pi, s, pjs, cj] == pj, pi == pis, pj == pjs}, {{pi, 0.3}, {pj, 0.4}, {pis, 0.35}, {pjs, 0.45}}] pjsol[s_Real?Positive] := pj /. FindRoot[{pbr[pj, s, pis, ci] == pi, pbr[pi, s, pjs, cj] == pj, pi == pis, pj == pjs}, {{pi, 0.3}, {pj, 0.4}, {pis, 0.35}, {pjs, 0.45}}] pisol[0.25] Plot[{pisol[s], pjsol[s]}, {s, 0., 0.5}](*Does not solve in 45 min. *)  I have tried to modify this to use derivatives instead of NArgMax: (*Using derivatives instead of NArgMax*) Clear[d, s, vi, vj, vlb, \ mui, muj, pi, pj, profit, ci, cj, pfoc, pisol, pjsol, pis, pjs]$Assumptions =
0 < {s, pi, pj, pis, pjs} < 1/2 && 0 < {mui, muj, vlb, vi, vj} < 1 &&
0 <= {ci, cj} < 1;
vlb[pi_, pj_] =
1 + pi - pj -
Sqrt Sqrt[s]; ci = 0.; cj = 0.05; mui = 0.5; muj = 0.5;
d[pi_, pj_, s_, pis_] =
mui*(1 - Max[pi, vlb[pi, pj]]) +
mui*Integrate[vi - pi + pj, {vi, Max[pi, 0], Max[pi, vlb[pi, pj]]},
Assumptions -> 0 < {s, pi, pj, pis, pjs} < 1/2 && 0 < {vlb} < 1] +
muj*Integrate[(1 - Max[pi, vj - pj + pi]), {vj, 0, vlb[pj, pis]},
Assumptions -> 0 < {s, pi, pj, pis, pjs} < 1/2 && 0 < {vlb} < 1];
d[0.2, 0.3, 0.1, 0.25]
profit[pi_, pj_, s_, pis_, c_] = (pi - c)*d[pi, pj, s, pis];
pfoc[pi_Real?Positive, pj_Real?Positive, s_Real?Positive,
pis_Real?Positive, c_Real?Positive] =
D[profit[pi, pj, s, pis, c], pi]
pisol[s_Real?Positive] :=
pi /. FindRoot[{pfoc[pi, pj, s, pis, ci] == 0,
pfoc[pj, pi, s, pjs, cj] == 0, pi == pis,
pj == pjs}, {{pi, 0.3}, {pj, 0.4}, {pis, 0.35}, {pjs, 0.45}}]
pjsol[s_Real?Positive] :=
pj /. FindRoot[{pfoc[pi, pj, s, pis, ci] == 0,
pfoc[pj, pi, s, pjs, cj] == 0, pi == pis,
pj == pjs}, {{pi, 0.3}, {pj, 0.4}, {pis, 0.35}, {pjs, 0.45}}]
pisol[0.25]
Plot[{pisol[s], pjsol[s]}, {s, 0., 0.5}]


Integrate was very slow on the Max of linear functions. The result should be piecewise quadratic and can be done by hand. NIntegrate gives errors about non-numerical boundaries even when replacing the symbolic boundary by a number:

Clear[d, s, vi, vj, vlb, mui, muj, pi, pj, profit, ci, cj, pfoc, \
pisol, pjsol, pis, pjs]
$Assumptions = 0 < {s, pi, pj, pis, pjs} < 1/2 && 0 < {mui, muj, vlb, vi, vj} < 1 && 0 <= {ci, cj} < 1; vlb[pi_, pj_] = 1 + pi - pj - Sqrt Sqrt[s]; ci = 0.; cj = 0.05; mui = 0.5; muj = 0.5; d[pi_Real?Positive, pj_Real?Positive, s_Real?Positive, pis_Real?Positive] = mui*(1 - Max[pi, vlb[pi, pj]]) + mui*NIntegrate[ vi - pi + pj, {vi, Max[pi, 0], Max[pi, vlb[pi, pj]]}] + muj*NIntegrate[(1 - Max[pi, vj - pj + pi]), {vj, 0, vlb[pj, pis]}]; d[0.2, 0.3, 0.1, 0.25]  yields vi = Max[0.,pi] is not a valid limit of integration. I have encountered similar problems many times before in Mathematica when trying to find fixed points. I gave up each time after many days of trying, so I do not anticipate a solution now. As I see it, there is a conflict between Mathematica constraints: Integrate does not work (for no good reason) or takes hours, but NIntegrate requires numeric arguments. For all parameters, numeric values would be fine, but for the candidate fixed point (pis in the above code), the maximization or derivative-taking has to happen before setting pis==pi. One way may be to solve using all possible guesses of pis before pis==pi, maybe by FunctionInterpolation, but this again is too slow. • Why are you using := and not =? Clear[x, y, xsola, ysola]; {xsola, ysola} = {x, y} /. FindRoot[{D[(x - 1)^2] == 0, D[(y - 2)^2] == 0}, {{x, 0.3}, {y, 0.4}}] – J. M.'s technical difficulties Feb 2 at 14:35 • @J.M. is in limbo I am trying to use this in a bigger piece of code where = creates errors of non-numeric values: FindRoot: The function value {<<1>>} is not a list of numbers with dimensions {2} at {pi,pj} = {0.3,0.4} I can expand the MWE to illustrate if my explanation is unclear. Basically, replace (x-1) by x-a and try to plot on {a,0,1}. – Sander Heinsalu Feb 2 at 14:46 • That looks like a case where FindRoot[] is not able to come up with a solution, so of course a replacement rule does not get produced. – J. M.'s technical difficulties Feb 2 at 14:54 • Could you clarify? FindRoot does find solutions in my MWE and the real code. Just assigning these solutions seems problematic. I can assign these element by element and plot just fine, but it seems inelegant. I can live with not assigning these as vectors - one of Mathematica's little mysteries why vector assignment does not work. – Sander Heinsalu Feb 2 at 15:02 • As I already said, your use of := is precisely the problem. Recall that := does not immediately evaluate its RHS, so of course there is a mismatch between {xsola, ysola} and ReplaceAll[{x, y}, FindRoot[(* stuff *)]]. Otherwise, the only way an assignment with = (instead of :=) will fail is if FindRoot[] itself cannot find a solution. Perhaps this is a time to reconsider the design of your actual code. – J. M.'s technical difficulties Feb 2 at 15:11 ## 2 Answers The offending line with the SetDelayed::shape error is: {xsola, ysola} := {x, y} /. FindRoot[{D[(x - 1)^2] == 0, D[(y - 2)^2] == 0}, {{x, 0.3}, {y, 0.4}}]  The FullForm of the offending line looks like SetDelay[List[..], ReplaceAll[..]]. The two arguments to SetDelayed[] should be lists of the same shape, if the first argument is a list. In this case, the second argument is not a list at all. When the first argument is a list, SetDelayed makes a series of definitions, one for each element in the list. The value for each comes from the corresponding element in the second argument. (It does this recursively down nest lists.) With Set instead of SetDelayed as @J.M. recommended, the ReplaceAll argument will evaluate to a list. Given this tidbit: I am trying to use this in a bigger piece of code where = creates errors of non-numeric values. It sounds like the problem is to delay the call of FindRoot until certain conditions are met. It might be better to refactor the code so that those conditions can be passed in functions calls. Not knowing more, something like the following might be a quick workaround: ClearAll[getSola]; getSola[] := {xsola, ysola} = {x, y} /. FindRoot[{D[(x - 1)^2] == 0, D[(y - 2)^2] == 0}, {{x, 0.3}, {y, 0.4}}]; getSola[] {xsol, ysol, xsola, ysola} (* {1., 2.} {1., 2., 1., 2.} *)  Also, one could hold up evaluation as follows depending on a condition, or give an error if the condition was not satisfied: ClearAll[getSola]; SetAttributes[getSola, HoldAll]; getSola[vars_, cond_] /; TrueQ@cond := vars = {x, y} /. FindRoot[{D[(x - a)^2] == 0, D[(y - b)^2] == 0}, {{x, 0.3}, {y, 0.4}}]; getSola[{xsola, ysola}, VectorQ[{a, b}, NumericQ]]  If a and b are numeric, FindRoot is called and the variables are set; otherwise, nothing is done. getSol[] is similar to ConditionalExpression[], except that ConditionalExpression[] evaluates to Undefined when the condition is False and getSola[] remains unevaluated. Probably the first alternative is better. At least it fits better with how I imagine the bigger code is organized. But I could easily be wrong. • I didn't want to guess at what monster the OP actually has, but your proposal is surely reasonable. – J. M.'s technical difficulties Feb 2 at 17:54 Long comment about an important side-issue in one of the examples: In one case, the assumptions are defined wrong: a < {x, y, z} < b does not translate into three inequalities. One can refactor with VectorLess or Thread: Clear[d, s, vi, vj, vlb, mui, muj, pi, pj, profit, ci, cj, pfoc, \ pisol, pjsol, pis, pjs] ($Assumptions =
Flatten@{Thread[0 < {s, pi, pj, pis, pjs} < 1/2],
Thread[0 < {mui, muj, vlb, vi, vj} < 1],
Thread[0 <= {ci, cj} < 1]};
vlb[pi_, pj_] = 1 + pi - pj - Sqrt Sqrt[s]; ci = 0.; cj = 0.05;
mui = 0.5; muj = 0.5;
d[pi_Real?Positive, pj_Real?Positive, s_Real?Positive,
pis_Real?Positive] =
mui*(1 - Max[pi, vlb[pi, pj]]) +
mui*Integrate[
vi - pi + pj, {vi, Max[pi, 0], Max[pi, vlb[pi, pj]]}] +
muj*Integrate[(1 - Max[pi, vj - pj + pi]), {vj, 0,
vlb[pj, pis]}];
d[0.2, 0.3, 0.1, 0.25]) // AbsoluteTiming

(* {0.162937, 0.545695}  *)


Note that exact (symbolic) solvers sometimes run into trouble with inexact floating-point parameters like ci = 0.1. Here it does not, so we're lucky. However, if it did, then one could do the following (and it's faster by a factor of 1/3000 if you have to call d[] multiple times): If you define d[] before you define the parameters, then the integration will be done with exact symbols just as quickly; and d[] returns the same result. Then you don't have to worry about round-off error causing Integrate to choke.

• Thank you for the assumption threading - I did not know it and was unaware of the issue, so could not even ask the question. The floating point warning messages with exact solvers I have seen before, but in this case did not notice - thanks for correcting that as well! – Sander Heinsalu Feb 3 at 2:38
• @SanderHeinsalu You're welcome! – Michael E2 Feb 3 at 2:41