How to solve a system of non-linear equations with known constants and unknown constants?

Question

I review these questions Q1 and Q2.

I wrote the following codes:

DEL = (1/2258332876800)*(80080*Subscript[c, 1]^6*Subscript[d, 1]^2 - 112*Subscript[c, 1]^5*Subscript[d, 1]*
     (-147900 + 46189*Subscript[d, 1] + 2145*Subscript[d, 1]^2 + 9724*Subscript[d, 2]) + 
    4*Subscript[c, 1]^4*(2246244*Subscript[d, 1]^3 + 45045*Subscript[d, 1]^4 - 1547*Subscript[d, 1]^2*
       (-10995 + 352*Subscript[c, 2] - 528*Subscript[d, 2]) + 23800*Subscript[d, 1]*(-6617 + 693*Subscript[d, 2]) + 
      272*(1147861 - 111930*Subscript[d, 2] + 3465*Subscript[d, 2]^2)) - 
    68*Subscript[c, 1]^3*(27027*Subscript[d, 1]^4 - 840*Subscript[d, 1]^2*(-1261 + 1452*Subscript[c, 2] - 1881*Subscript[d, 2]) + 
      Subscript[d, 1]^3*(765135 - 48048*Subscript[c, 2] + 36036*Subscript[d, 2]) + 
      9984*(18513 - 6160*Subscript[d, 2] + 315*Subscript[d, 2]^2) - 16*Subscript[d, 1]*(2024737 - 876330*Subscript[d, 2] - 
        10395*Subscript[d, 2]^2 + 280*Subscript[c, 2]*(-884 + 99*Subscript[d, 2]))) + 
    1400256*(336*(27*Subscript[d, 1]^2 + 24*Subscript[d, 1]*(5 + 3*Subscript[d, 2]) + 
        16*(15 + 10*Subscript[d, 2] + 3*Subscript[d, 2]^2)) - 24*Subscript[c, 2]*(135*Subscript[d, 1]^2 + 
        24*Subscript[d, 1]*(-14 + 15*Subscript[d, 2]) + 16*(-175 - 28*Subscript[d, 2] + 15*Subscript[d, 2]^2)) + 
      Subscript[c, 2]^2*(315*Subscript[d, 1]^2 + 24*Subscript[d, 1]*(-81 + 35*Subscript[d, 2]) + 
        16*(903 - 162*Subscript[d, 2] + 35*Subscript[d, 2]^2))) + 17*Subscript[c, 1]^2*(887040*Subscript[c, 2]^2*Subscript[d, 1]^2 + 
      280665*Subscript[d, 1]^4 + 68040*Subscript[d, 1]^3*(117 + 11*Subscript[d, 2]) + 
      6150144*(603 - 234*Subscript[d, 2] + 35*Subscript[d, 2]^2) + 139776*Subscript[d, 1]*(-693 + 1430*Subscript[d, 2] + 
        135*Subscript[d, 2]^2) + 16*Subscript[d, 1]^2*(6961669 + 1547910*Subscript[d, 2] + 31185*Subscript[d, 2]^2) - 
      96*Subscript[c, 2]*(20790*Subscript[d, 1]^3 + 35*Subscript[d, 1]^2*(7709 + 792*Subscript[d, 2]) + 
        520*Subscript[d, 1]*(-3784 + 1197*Subscript[d, 2]) + 208*(18513 - 6160*Subscript[d, 2] + 315*Subscript[d, 2]^2))) - 
    10608*Subscript[c, 1]*(80*Subscript[c, 2]^2*Subscript[d, 1]*(-484 + 189*Subscript[d, 1] + 252*Subscript[d, 2]) + 
      132*(315*Subscript[d, 1]^3 + 24*Subscript[d, 1]^2*(57 + 35*Subscript[d, 2]) + 
        768*(-175 - 28*Subscript[d, 2] + 15*Subscript[d, 2]^2) + 16*Subscript[d, 1]*(-1869 + 654*Subscript[d, 2] + 
          35*Subscript[d, 2]^2)) - Subscript[c, 2]*(8505*Subscript[d, 1]^3 + 840*Subscript[d, 1]^2*(110 + 27*Subscript[d, 2]) + 
        8448*(903 - 162*Subscript[d, 2] + 35*Subscript[d, 2]^2) + 16*Subscript[d, 1]*(-14949 + 21560*Subscript[d, 2] + 
          945*Subscript[d, 2]^2))))

I want to solve the following system equations which Subscript[c, 1],Subscript[c, 2 are unknown variables and Subscript[d, 1],Subscript[d, 2 are known variables:

E1 = D[DEL, {Subscript[c, 1], 1}]; 
E2 = D[DEL, {Subscript[c, 2], 1}]; 
Solve[E1 == 0 && E2 == 0, {Subscript[c, 1], Subscript[c, 2]}, Reals]

Any suggestions?

Answer 1

I wanted to inform the system that Subscript[d,1] and Subscript[d,2] were constants. I was unable to do that because, for example, Subscript[d,1] is not a symbol (one often runs into problems trying to use subscripted symbols in Mathematica).

So I modified your definition of DEL and replace all subscripted variables with a non-subscripted version.

Using your definition of DEL (not repeated here).

DEL = DEL /. Subscript[a_, i_] :> ToExpression[ToString[a] ~~ ToString[i]]

SetAttributes[d1, Constant]
SetAttributes[d2, Constant]

Then proceed as in your query using Daniel Lichtblau's tip about not constraining the result to reals.

E1 = D[DEL, {c1, 1}];
E2 = D[DEL, {c2, 1}];

Now Solve produces a result (it is verbose, but at least I get a result)

Solve[E1 == 0 && E2 == 0, {c1, c2}]

In the result on my system there were nine solutions.

I copied the first solution (again, very verbose) and created the functions

c1[d1_, d2_] := N["copy of solution for c1"]
c2[d1_, d2_] := N["copy of solution for c2"]

Then I could, for example, plot c1[x,y] in a 3D plot

Plot3D[c1[x, y], {x, -10, 10}, {y, -10, 10}]

The second function was slow so I made a table of values and used ListPlot3D.

table = Flatten[Table[{x, y, c2[x, y]}, {x, -10, 10}, {y, -10, 10}], 
   1];

ListPlot3D[table]

Answer 2

You can try to find numerical solution using FindRoot by fixing d1 and d2.

Subscript[d, 2] = 1; Subscript[d, 1] = 1;
Plot3D[{E1, E2}, {Subscript[c, 1], -10, 10}, {Subscript[c, 2], -10,10}]

ContourPlot[{E1, E2}, {Subscript[c, 1], -5, 5}, {Subscript[c, 2], -12,12}]

Reap[FindRoot[{E1 == 0 && E2 == 0}, {{Subscript[c, 1], 4}, {Subscript[c, 2], 4}}, 
  StepMonitor :> Sow[{Subscript[c, 1], Subscript[c, 2]}]]]

{{Subscript[c, 1] -> 4.01292, Subscript[c, 2] -> -9.97329}, {{{3.32569, -9.05764}, {4.33577, -10.5018}, \ {4.05056, -10.032}, {4.01351, -9.97422}, {4.01292, -9.97329}, \ {4.01292, -9.97329}, {4.01292, -9.97329}}}}