0
$\begingroup$

The following code symbolically computes a system of equations and then applies FindRoot to it. When I run code the message 'Encountered a singular Jacobian...' comes up, for all initial points that I have tried.

k = 0.5;
fA = Max[0, Min[1, (1/h)*(Log[(pAC - k)/(pA - k)] + pAC/(pAC - k) - pA*(pA - k))]]; 
fB = Max[0, Min[1, (1/h)*(Log[(pBC - k)/(pB - k)] + pBC/(pBC - k) - pB*(pB - k))]]; 
zAB = Max[0, Min[1, 1/2 + (1/(2*g))*(Log[(pB - k)/(pA - k)] + pB/(pB - k) - pA/(pA - k))]]; 
zC = Max[0, Min[1, 1/2 + (1/(2*g))*(Log[(pBC - k)/(pAC - k)] + pBC/(pBC - k) - pAC/(pAC - k))]]; 
rAA = fA*zAB - (1/2)*Max[0, zAB - zC]*Max[0, fA - fB]; 
rBB = fB*(1 - zAB) - (1/2)*Max[0, (1 - zAB) - (1 - zC)]*Max[0, fB - fA]; 
rAC = (1 - fA)*zC - (1/2)*Max[0, zC - zAB]*Max[0, (1 - fA) - (1 - fB)]; 
rBC = (1 - fB)*(1 - zC) - (1/2)*Max[0, (1 - zC) - (1 - zAB)]*Max[0, (1 - fB) - (1 - fA)]; 
disutilityAA = fA*zAB*(h*(fA/2) + g*(zAB/2)) - (1/2)*Max[0, zAB - zC]*Max[0, fA - fB]*(h*(fA - (1/3)*(fA - fB)) + g*(zAB - (1/3)*(zAB - zC))); 
disutilityBB = fB*(1 - zAB)*(h*(fB/2) + g*((1 - zAB)/2)) - (1/2)*Max[0, (1 - zAB) - (1 - zC)]*Max[0, fB - fA]*
     (h*(fB - (1/3)*(fB - fA)) + g*((1 - zAB) - (1/3)*((1 - zAB) - (1 - zC)))); 
disutilityAC = (1 - fA)*zC*(h*((1 - fA)/2) + g*(zC/2)) - (1/2)*Max[0, zC - zAB]*Max[0, (1 - fA) - (1 - fB)]*
     (h*((1 - fA) - (1/3)*((1 - fA) - (1 - fB)) + g*(zC - (1/3)*(zC - zAB)))); 
disutilityBC = (1 - fB)*(1 - zC)*(h*((1 - fB)/2) + g*((1 - zC)/2)) - (1/2)*Max[0, (1 - zC) - (1 - zAB)]*(h*((1 - fB) - (1/3)*((1 - fB) - (1 - fA))) + 
      g*((1 - zC) - (1/3)*((1 - zC) - (1 - zAB))));
everythingintermsofplayersactions = {xA -> 1/(pA - k), xAC -> 1/(pAC - k), xB -> 1/(pB - k), xBC -> 1/(pBC - k)};
\[Pi]A = (pA/(1 + tA) - c)*xA*(mA + rAA*mC) + (pAC/(1 + tC) - c)*xAC*rAC*mC //. everythingintermsofplayersactions; 
\[Pi]B = (pB/(1 + tB) - c)*xB*(mB + rBB*mC) + (pBC/(1 + tC) - c)*xBC*rBC*mC //. everythingintermsofplayersactions; 
focA = D[\[Pi]A, pA] == 0; 
focAC = D[\[Pi]A, pAC] == 0; 
focB = D[\[Pi]B, pB] == 0; 
focBC = D[\[Pi]B, pBC] == 0; 
dependenceontaxes = {pA -> pA[tA, tB, tC], pB -> pB[tA, tB, tC]}; 
focsp = {focA, focAC, focB, focBC} /. dependenceontaxes; 
t = {tA, tB, tC}; 
dfocspbydt = D[focsp, {t}]; 
\[Pi]Aasafunctionoft = \[Pi]A /. dependenceontaxes; 
\[Pi]Basafunctionoft = \[Pi]B /. dependenceontaxes; 
sA = \[Pi]A + (tA/(1 + tA))*pA*xA*(mA + mC*rAA) + mA*(xA*k - xA*pA) //. everythingintermsofplayersactions /. dependenceontaxes; 
sB = \[Pi]B + (tB/(1 + tB))*pB*xB*(mB + mC*rBB) + mB*(yB*k - xB*pB) //. everythingintermsofplayersactions /. dependenceontaxes; 
  sC = (tC/(1 + tC))*mC*(rBC*xBC*pBC + rAC*xAC*pAC) + mC*(k*(rAA*xA + rBB*xB + rAC*xAC + rBC*xBC) - rAA*xA*pA - rBB*xB*pB - rAC*xAC*pAC - rBC*xBC*pBC) - 
      disutilityAA - disutilityAC - disutilityBB - disutilityBC //. everythingintermsofplayersactions /. dependenceontaxes; 
focsA = D[sA, tA] == 0; 
focsB = D[sB, tB] == 0; 
focsC = D[sC, tC] == 0; 
allconditions = Join[Flatten[focsp], Flatten[dfocspbydt], {focsA}, {focsB}, {focsC}]; 
allconditionswithoutfunctions = 
   allconditions /. pAC[tA, tB, tC] -> pAC /. pA[tA, tB, tC] -> pA /. pB[tA, tB, tC] -> pB /. pBC[tA, tB, tC] -> pBC /. Derivative[1, 0, 0][pA][tA, tB, tC] -> 
                dpAdtA /. Derivative[0, 1, 0][pA][tA, tB, tC] -> dpAdtB /. Derivative[0, 0, 1][pA][tA, tB, tC] -> dpAdtC /. Derivative[1, 0, 0][pAC][tA, tB, tC] -> 
             dpACdtA /. Derivative[0, 1, 0][pAC][tA, tB, tC] -> dpACdtB /. Derivative[0, 0, 1][pAC][tA, tB, tC] -> dpACdtC /. 
         Derivative[1, 0, 0][pB][tA, tB, tC] -> dpBdtA /. Derivative[0, 1, 0][pB][tA, tB, tC] -> dpBdtB /. Derivative[0, 0, 1][pB][tA, tB, tC] -> dpBdtC /. 
      Derivative[1, 0, 0][pBC][tA, tB, tC] -> dpBCdtA /. Derivative[0, 1, 0][pBC][tA, tB, tC] -> dpBCdtB /. Derivative[0, 0, 1][pBC][tA, tB, tC] -> dpBCdtC; 
parameters = {\[Rho] -> 1/2, c -> 1, mA -> 1/3, mB -> 1/3, mC -> 1/3, g -> 0.15, h -> 0.2}; 
allconditionswithoutfunctionswithnumbers = allconditionswithoutfunctions /. parameters; 
FindRoot[{allconditionswithoutfunctionswithnumbers}, {{tA, 0.43}, {tB, 0.44}, {tC, 0.88}, {pAC, 6.1}, {pA, 4.}, {pBC, 7.1}, {pB, 3.2}, {dpAdtA, 1.1}, {dpAdtB, 1.2}, 
   {dpAdtC, 1.05}, {dpACdtA, 1.06}, {dpACdtB, 1.03}, {dpACdtC, 1.02}, {dpBdtA, 1.04}, {dpBdtB, 1.03}, {dpBdtC, 1.01}, {dpBCdtA, 1.02}, {dpBCdtB, 1.04}, 
   {dpBCdtC, 1.03}}]

The system consists of 19 equations and has 19 variables. Is it possible that this is due to a problem with FindRoot? (From the context of the problem it seems that there should be a solution and that the solutions should be locally unique.)

When I change the value of k to k=0 at the beginning, then I instead always get the error message 'The line search decreased the step size to within tolerance specified \ by AccuracyGoal and PrecisionGoal but was unable to find a sufficient \ decrease in the merit function. You may need more than \ MachinePrecision digits of working precision to meet these tolerances.'

$\endgroup$
  • $\begingroup$ How does Sard's theorem collide with the fact that the mapping $f \colon \mathbb{R}^n \to \mathbb{R}^n$, $f(x)=0$ has noninvertible differential? Not at all. $\endgroup$ – Henrik Schumacher Mar 3 '18 at 14:27
  • 1
    $\begingroup$ Sard's theorem only says that the set of all critical values is a null set. By definition, critical values are within the image of the mapping. $\endgroup$ – Henrik Schumacher Mar 3 '18 at 14:29
  • $\begingroup$ Sorry, yes, I was mistaken. I have edited the question to correct it. $\endgroup$ – Lennart Mar 3 '18 at 14:41
  • 1
    $\begingroup$ That does not suffice for Newton's algorithm (which is internally used in FindRoot). In patricular, I think the line search relies on the fact that the expression is comntinuously differentiable. There are algorithms that can deal with, e.g., semismooth expressions but I doubt that Mathematica employs such an algorithm; it just does not know how to compute a (generalized) derivative. However, you can help Mathematica by supplying a costumized Jacobian via the option Jacobian... $\endgroup$ – Henrik Schumacher Mar 3 '18 at 15:01
  • 2
    $\begingroup$ The use of Min and Max means there may be "flat" regions, and this will manifest as regions where the Jacobian is rank deficient. $\endgroup$ – Daniel Lichtblau Mar 3 '18 at 16:23
5
$\begingroup$

In the meantime, I did some experiments. I repaced all occurences of Max[0,x] by Ramp[x] and Min[1,x] by 1 - Ramp[1 - x] since these can be easily differentiated by Mathematica.

Moreover, I computed the derivative of your expression

Off[Part::partd];
vars = {tA, tB, tC, pAC, pA, pBC, pB, dpAdtA, dpAdtB, dpAdtC, dpACdtA,
    dpACdtB, dpACdtC, dpBdtA, dpBdtB, dpBdtC, dpBCdtA, dpBCdtB, 
   dpBCdtC};
XX = Table[X[[i]], {i, 1, Length[vars]}];
expr = allconditionswithoutfunctionswithnumbers[[All, 1]];
F = X \[Function] Evaluate[N[expr /. Thread[vars -> XX]]];
DF = X \[Function] Evaluate[N[D[F[XX], {XX, 1}]]];

Let's calculate the derivative at the initial point and check its rank:

X0 = {0.43, 0.44, 0.88, 6.1, 4., 7.1, 3.2, 1.1, 1.2, 1.05, 1.06, 1.03,
    1.02, 1.04, 1.03, 1.01, 1.02, 1.04, 1.03};
Dimensions[DF[X0]]
MatrixRank[DF[X0]]

{19, 19}

12

Nope, it's not invertible, even with lots of good will.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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