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.
:=
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$=
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)
byx-a
and try to plot on{a,0,1}
. $\endgroup$FindRoot[]
is not able to come up with a solution, so of course a replacement rule does not get produced. $\endgroup$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$:=
is precisely the problem. Recall that:=
does not immediately evaluate its RHS, so of course there is a mismatch between{xsola, ysola}
andReplaceAll[{x, y}, FindRoot[(* stuff *)]]
. Otherwise, the only way an assignment with=
(instead of:=
) will fail is ifFindRoot[]
itself cannot find a solution. Perhaps this is a time to reconsider the design of your actual code. $\endgroup$