0
$\begingroup$

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[2] 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[2] 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[2] 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.

$\endgroup$
8
  • 1
    $\begingroup$ 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}}] $\endgroup$ Feb 2, 2020 at 14:35
  • $\begingroup$ @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}. $\endgroup$ Feb 2, 2020 at 14:46
  • $\begingroup$ 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. $\endgroup$ Feb 2, 2020 at 14:54
  • $\begingroup$ 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. $\endgroup$ Feb 2, 2020 at 15:02
  • $\begingroup$ 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. $\endgroup$ Feb 2, 2020 at 15:11

2 Answers 2

5
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ I didn't want to guess at what monster the OP actually has, but your proposal is surely reasonable. $\endgroup$ Feb 2, 2020 at 17:54
2
$\begingroup$

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[2] 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.

$\endgroup$
2
  • $\begingroup$ 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! $\endgroup$ Feb 3, 2020 at 2:38
  • $\begingroup$ @SanderHeinsalu You're welcome! $\endgroup$
    – Michael E2
    Feb 3, 2020 at 2:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.