Skip to main content

equation solving - Solve answer optimization


I have an equation system (here is the simplified version for illustration)


    eq = Thread[
Equal[{a[

0]^2 (a[47]^2 (2 - 2 u[2] - 2 u[5] + 6 u[7] - 2 u[9] +
6 u[10] - 2 u[11] + 6 u[13] - 2 u[15] - 2 u[17]) +
a[49]^2 (2 - 2 u[2] - 2 u[5] + 6 u[7] - 2 u[9] + 6 u[10] -
2 u[11] + 6 u[13] - 2 u[15] - 2 u[17]) +
a[52]^2 (2 - 2 u[2] - 2 u[5] + 6 u[7] - 2 u[9] + 6 u[10] -
2 u[11] + 6 u[13] - 2 u[15] - 2 u[17])) +
a[47]^2 (a[
49]^2 (2 + 2 u[1] - 2 u[2] + 2 u[3] - 6 u[4] + 6 u[5] -
6 u[6] - 2 u[7] + 2 u[8] - 2 u[9] + 6 u[10] - 2 u[11] +
2 u[12] - 2 u[13] - 6 u[14] + 6 u[15] + 2 u[16] -

2 u[17] + 2 u[18] - 6 u[19]) +
a[52]^2 (2 + 2 u[1] - 2 u[2] + 2 u[3] - 6 u[4] + 6 u[5] -
6 u[6] - 2 u[7] + 2 u[8] - 2 u[9] + 6 u[10] - 2 u[11] +
2 u[12] - 2 u[13] - 6 u[14] + 6 u[15] + 2 u[16] -
2 u[17] + 2 u[18] - 6 u[19])) +
a[0] a[47] a[49] a[
52] (8 - 12 u[1] + 12 u[3] + 12 u[4] - 12 u[6] + 12 u[8] -
24 u[9] + 24 u[10] + 24 u[11] - 12 u[12] + 12 u[14] -
12 u[16] + 12 u[18] - 12 u[19]) +
a[49]^2 a[

52]^2 (2 + 2 u[1] - 2 u[2] + 2 u[3] - 6 u[4] + 6 u[5] -
6 u[6] - 2 u[7] + 2 u[8] - 2 u[9] + 6 u[10] - 2 u[11] +
2 u[12] - 2 u[13] - 6 u[14] + 6 u[15] + 2 u[16] - 2 u[17] +
2 u[18] - 6 u[19]) +
a[47]^4 (-1 - u[1] + u[2] - u[3] - u[4] + u[5] - u[6] + u[7] -
u[8] + u[9] + u[10] + u[11] - u[12] + u[13] - u[14] + u[15] -
u[16] + u[17] - u[18] - u[19]) +
a[49]^4 (-1 - u[1] + u[2] - u[3] - u[4] + u[5] - u[6] + u[7] -
u[8] + u[9] + u[10] + u[11] - u[12] + u[13] - u[14] + u[15] -
u[16] + u[17] - u[18] - u[19]) +

a[52]^4 (-1 - u[1] + u[2] - u[3] - u[4] + u[5] - u[6] + u[7] -
u[8] + u[9] + u[10] + u[11] - u[12] + u[13] - u[14] + u[15] -
u[16] + u[17] - u[18] - u[19]) +
a[0]^4 (-1 + u[1] + u[2] + u[3] + u[4] + u[5] + u[6] + u[7] +
u[8] + u[9] + u[10] + u[11] + u[12] + u[13] + u[14] + u[15] +
u[16] + u[17] + u[18] + u[19]),
a[0]^2 a[49] a[
52] (-8 u[1] + 4 u[2] + 8 u[3] + 12 u[4] - 4 u[6] - 4 u[7] -
4 u[9] + 8 u[11] - 4 u[12] - 8 u[13] + 4 u[15] + 4 u[16] +
4 u[18]) +

a[0]^3 a[
47] (2 u[1] - 2 u[3] - 2 u[4] + 2 u[6] + 4 u[7] + 2 u[8] +
4 u[10] + 2 u[12] + 4 u[13] + 2 u[14] - 2 u[16] + 2 u[18] -
2 u[19]) +
a[49]^3 a[
52] (4 u[4] - 4 u[5] + 4 u[6] - 4 u[10] + 4 u[14] - 4 u[15] +
4 u[19]) +
a[49] a[52]^3 (4 u[4] - 4 u[5] + 4 u[6] - 4 u[10] + 4 u[14] -
4 u[15] + 4 u[19]) +
a[47]^2 a[49] a[

52] (4 u[1] - 4 u[2] + 4 u[3] - 4 u[7] + 4 u[8] - 4 u[9] +
8 u[11] - 8 u[12] + 8 u[13] + 12 u[14] - 12 u[15] -
8 u[16] + 8 u[17] - 8 u[18] + 12 u[19]) +
a[0] (a[47]^3 (-2 u[1] + 2 u[3] + 2 u[4] - 2 u[6] + 4 u[7] -
2 u[8] + 4 u[10] - 2 u[12] + 4 u[13] - 2 u[14] +
2 u[16] - 2 u[18] + 2 u[19]) +
a[47] (a[
49]^2 (6 u[1] - 4 u[2] + 2 u[3] - 6 u[4] + 4 u[5] -
2 u[6] - 2 u[8] + 4 u[9] - 8 u[11] + 2 u[12] +
4 u[13] - 6 u[14] + 6 u[16] - 6 u[18] + 6 u[19]) +

a[52]^2 (6 u[1] - 4 u[2] + 2 u[3] - 6 u[4] + 4 u[5] -
2 u[6] - 2 u[8] + 4 u[9] - 8 u[11] + 2 u[12] + 4 u[13]
- 6 u[14] + 6 u[16] - 6 u[18] + 6 u[19]))),
a[0]^2 a[47] a[
52] (-8 u[1] + 4 u[2] + 8 u[3] + 12 u[4] - 4 u[6] - 4 u[7] -
4 u[9] + 8 u[11] - 4 u[12] - 8 u[13] + 4 u[15] + 4 u[16] +
4 u[18]) +
a[0]^3 a[
49] (2 u[1] - 2 u[3] - 2 u[4] + 2 u[6] + 4 u[7] + 2 u[8] +
4 u[10] + 2 u[12] + 4 u[13] + 2 u[14] - 2 u[16] + 2 u[18] -

2 u[19]) +
a[47]^3 a[
52] (4 u[4] - 4 u[5] + 4 u[6] - 4 u[10] + 4 u[14] - 4 u[15] +
4 u[19]) +
a[47] a[52]^3 (4 u[4] - 4 u[5] + 4 u[6] - 4 u[10] + 4 u[14] -
4 u[15] + 4 u[19]) +
a[47] a[49]^2 a[
52] (4 u[1] - 4 u[2] + 4 u[3] - 4 u[7] + 4 u[8] - 4 u[9] +
8 u[11] - 8 u[12] + 8 u[13] + 12 u[14] - 12 u[15] -
8 u[16] + 8 u[17] - 8 u[18] + 12 u[19]) +

a[0] (a[49]^3 (-2 u[1] + 2 u[3] + 2 u[4] - 2 u[6] + 4 u[7] -
2 u[8] + 4 u[10] - 2 u[12] + 4 u[13] - 2 u[14] +
2 u[16] - 2 u[18] + 2 u[19]) +
a[49] (a[
47]^2 (6 u[1] - 4 u[2] + 2 u[3] - 6 u[4] + 4 u[5] -
2 u[6] - 2 u[8] + 4 u[9] - 8 u[11] + 2 u[12] +
4 u[13] - 6 u[14] + 6 u[16] - 6 u[18] + 6 u[19]) +
a[52]^2 (6 u[1] - 4 u[2] + 2 u[3] - 6 u[4] + 4 u[5] -
2 u[6] - 2 u[8] + 4 u[9] - 8 u[11] + 2 u[12] +
4 u[13] - 6 u[14] + 6 u[16] - 6 u[18] + 6 u[19]))),

a[0]^2 a[47] a[
49] (-8 u[1] + 4 u[2] + 8 u[3] + 12 u[4] - 4 u[6] - 4 u[7] -
4 u[9] + 8 u[11] - 4 u[12] - 8 u[13] + 4 u[15] + 4 u[16] +
4 u[18]) +
a[0]^3 a[
52] (2 u[1] - 2 u[3] - 2 u[4] + 2 u[6] + 4 u[7] + 2 u[8] +
4 u[10] + 2 u[12] + 4 u[13] + 2 u[14] - 2 u[16] + 2 u[18] -
2 u[19]) +
a[47]^3 a[
49] (4 u[4] - 4 u[5] + 4 u[6] - 4 u[10] + 4 u[14] - 4 u[15] +

4 u[19]) +
a[0] (a[52]^3 (-2 u[1] + 2 u[3] + 2 u[4] - 2 u[6] + 4 u[7] -
2 u[8] + 4 u[10] - 2 u[12] + 4 u[13] - 2 u[14] +
2 u[16] - 2 u[18] + 2 u[19]) +
a[47]^2 a[
52] (6 u[1] - 4 u[2] + 2 u[3] - 6 u[4] + 4 u[5] - 2 u[6] -
2 u[8] + 4 u[9] - 8 u[11] + 2 u[12] + 4 u[13] - 6 u[14] +
6 u[16] - 6 u[18] + 6 u[19]) +
a[49]^2 a[
52] (6 u[1] - 4 u[2] + 2 u[3] - 6 u[4] + 4 u[5] - 2 u[6] -

2 u[8] + 4 u[9] - 8 u[11] + 2 u[12] + 4 u[13] - 6 u[14] +
6 u[16] - 6 u[18] + 6 u[19])) +
a[47] (a[
49]^3 (4 u[4] - 4 u[5] + 4 u[6] - 4 u[10] + 4 u[14] -
4 u[15] + 4 u[19]) +
a[49] a[52]^2 (4 u[1] - 4 u[2] + 4 u[3] - 4 u[7] + 4 u[8] -
4 u[9] + 8 u[11] - 8 u[12] + 8 u[13] + 12 u[14] -
12 u[15] - 8 u[16] + 8 u[17] - 8 u[18] + 12 u[19]))}, 0]];

I want to solve them for any parameters a[_]. For this simple equation the answer is



SolveAlways[eq, {a[0], a[47], a[49], a[52]}]

Out[2]= {{u[2] -> -u[4] - u[6] - u[8] - 2 u[12] + 2 u[13] - 2 u[15] - u[16] + u[17] - 3 u[18] + 2 u[19],


u[5] -> u[4] + u[6] - u[10] + u[14] - u[15] + u[19],


u[7] -> -u[10] - u[13],


u[9] -> 2/3 + u[8] + u[10] + u[12] - u[13] + u[15] - u[17] + 2 u[18] - 2 u[19],


u[11] -> 1/3 + u[12] - u[13] - u[14] + u[15] + u[16] - u[17] + u[18] - u[19],


u[1] -> -u[6] - u[8] - u[12] - u[14] - u[18],


u[3] -> -u[4] - u[16] - u[19]}}


So, I have 7 solution, from which only two of them have free constants 1/3 and 2/3. All other u[_], I can take zeros, so my simple answer is just u[9]->2/3 && u[11]->1/3



Unfortunately I need to solve much large system, which contains more that 7000 u[], and SolveAlways[ ] simply cannot do that. Fortunately, I can use linearity in u[] and avoid SolveAlways[] completely using simple variable replacement


eqToUspace = (eq /. {y_Plus?(FreeQ[#, a[_]] &) :> uSpace[y]});
(toReplace = Union[Cases[eqToUspace, _uSpace, Infinity]]);
newvarsList = Table[ur[i], {i, Length[toReplace]}];
uReplRules = Thread[Rule[toReplace, newvarsList]];

Which yields the system


simpleEq = eqToUspace /. uReplRules

Out[2]= {a[0]^2 (a[47]^2 ur[1] + a[49]^2 ur[1] + a[52]^2 ur[1]) + a[0] a[47] a[49] a[52] ur[3] + a[49]^2 a[52]^2 ur[4] + a[47]^2 (a[49]^2 ur[4] + a[52]^2 ur[4]) + a[47]^4 ur[6] + a[49]^4 ur[6] + a[52]^4 ur[6] + a[0]^4 ur[7] == 0, a[0]^2 a[49] a[52] ur[2] + a[0]^3 a[47] ur[5] + a[49]^3 a[52] ur[9] + a[49] a[52]^3 ur[9] + a[0] (a[47]^3 ur[8] + a[47] (a[49]^2 ur[10] + a[52]^2 ur[10])) + a[47]^2 a[49] a[52] ur[11] == 0, a[0]^2 a[47] a[52] ur[2] + a[0]^3 a[49] ur[5] + a[47]^3 a[52] ur[9] + a[47] a[52]^3 ur[9] + a[0] (a[49]^3 ur[8] + a[49] (a[47]^2 ur[10] + a[52]^2 ur[10])) + a[47] a[49]^2 a[52] ur[11] == 0, a[0]^2 a[47] a[49] ur[2] + a[0]^3 a[52] ur[5] + a[47]^3 a[49] ur[9] + a[0] (a[52]^3 ur[8] + a[47]^2 a[52] ur[10] + a[49]^2 a[52] ur[10]) + a[47] (a[49]^3 ur[9] + a[49] a[52]^2 ur[11]) == 0}



Solution of which (can be tested by SolveAlways[] as well) is just zeroes: {{ur[1] -> 0, ur[2] -> 0, ur[3] -> 0, ur[4] -> 0, ur[5] -> 0, ur[6] -> 0, ur[7] -> 0, ur[8] -> 0, ur[9] -> 0, ur[10] -> 0, ur[11] -> 0}}


After back substitution of old variables I can then solve the new system


unOptimized = 
Solve[Thread[Equal[toReplace /. {uSpace -> Identity}, 0]],
Cases[toReplace, _u, Infinity]]

During evaluation of In[]:= Solve::svars: Equations may not give solutions for all "solve" variables. >>


Out[167]= {{u[13] -> -u[7] - u[10], u[17] -> 1 - u[2] - u[5] - u[9] - u[11] - u[15], u[12] -> 1/3 - u[1] + u[3] + u[4] - u[6] - u[9] + u[10] + u[11], u[18] -> 1/3 + u[1] - u[2] - u[3] - 2 u[4] - u[7] - u[10] - u[11] - u[15] - u[16], u[8] -> -(2/3) - u[1] + u[2] - u[3] + u[4] - u[5] + u[6] + u[7] + u[9] - u[10], u[14] -> u[3] + u[5] - u[6] + u[10] + u[15] + u[16], u[19] -> -u[3] - u[4] - u[16]}}


It contains same number of solutions. The only problem is that it is unoptimal in the sense that instead of just 2 free constants I have 4. Of course this is related to which u[_] Solve decides to solve first.


My question is how to find optimal solution of the reformulated problem with smallest number of constant terms (i.e. which have no u[_]). I even cannot identify the type of problem here.




Answer



Here a way, may be not to reduce calculation time, but to reduce memory usage with large systems.


First define all possible parameter combinations of the a[i] and then find the corresponding equations for the u[i]. In order to be valid for all a[i], these equations have all to equal zero.


paracomb = (Outer[
Times, {a[0], a[47], a[49], a[52]}, {a[0], a[47], a[49],
a[52]}, {a[0], a[47], a[49], a[52]}, {a[0], a[47], a[49],
a[52]}] // Flatten // Union)

(* {a[0]^4, a[0]^3 a[47], a[0]^2 a[47]^2, a[0] a[47]^3, a[47]^4,
a[0]^3 a[49], a[0]^2 a[47] a[49], a[0] a[47]^2 a[49], a[47]^3 a[49],

a[0]^2 a[49]^2, a[0] a[47] a[49]^2, a[47]^2 a[49]^2, a[0] a[49]^3,
a[47] a[49]^3, a[49]^4, a[0]^3 a[52], a[0]^2 a[47] a[52],
a[0] a[47]^2 a[52], a[47]^3 a[52], a[0]^2 a[49] a[52],
a[0] a[47] a[49] a[52], a[47]^2 a[49] a[52], a[0] a[49]^2 a[52],
a[47] a[49]^2 a[52], a[49]^3 a[52], a[0]^2 a[52]^2,
a[0] a[47] a[52]^2, a[47]^2 a[52]^2, a[0] a[49] a[52]^2,
a[47] a[49] a[52]^2, a[49]^2 a[52]^2, a[0] a[52]^3, a[47] a[52]^3,
a[49] a[52]^3, a[52]^4} *)

coe = Coefficient[eq[[All, 1]], #, 1] & /@ paracomb;


coeu = Union[Flatten[DeleteCases[#, 0] & /@ coe] // Simplify,
SameTest -> (#1 === -#2 || #1 === #2 &)]

(* {-2 (-1 + u[2] + u[5] - 3 u[7] + u[9] - 3 u[10] + u[11] - 3 u[13] +
u[15] + u[17]), -4 (2 u[1] - u[2] - 2 u[3] - 3 u[4] + u[6] +
u[7] + u[9] - 2 u[11] + u[12] + 2 u[13] - u[15] - u[16] - u[18]),
2 (1 + u[1] - u[2] + u[3] - 3 u[4] + 3 u[5] - 3 u[6] - u[7] + u[8] -
u[9] + 3 u[10] - u[11] + u[12] - u[13] - 3 u[14] + 3 u[15] +
u[16] - u[17] + u[18] - 3 u[19]), -1 - u[1] + u[2] - u[3] - u[4] +

u[5] - u[6] + u[7] - u[8] + u[9] + u[10] + u[11] - u[12] + u[13] -
u[14] + u[15] - u[16] + u[17] - u[18] -
u[19], -2 (u[1] - u[3] - u[4] + u[6] - 2 u[7] + u[8] - 2 u[10] +
u[12] - 2 u[13] + u[14] - u[16] + u[18] - u[19]),
2 (u[1] - u[3] - u[4] + u[6] + 2 u[7] + u[8] + 2 u[10] + u[12] +
2 u[13] + u[14] - u[16] + u[18] - u[19]),
4 (u[4] - u[5] + u[6] - u[10] + u[14] - u[15] + u[19]), -1 + u[1] +
u[2] + u[3] + u[4] + u[5] + u[6] + u[7] + u[8] + u[9] + u[10] +
u[11] + u[12] + u[13] + u[14] + u[15] + u[16] + u[17] + u[18] +
u[19], -4 (-2 + 3 u[1] - 3 u[3] - 3 u[4] + 3 u[6] - 3 u[8] +

6 u[9] - 6 u[10] - 6 u[11] + 3 u[12] - 3 u[14] + 3 u[16] -
3 u[18] + 3 u[19]),
2 (3 u[1] - 2 u[2] + u[3] - 3 u[4] + 2 u[5] - u[6] - u[8] + 2 u[9] -
4 u[11] + u[12] + 2 u[13] - 3 u[14] + 3 u[16] - 3 u[18] +
3 u[19]),
4 (u[1] - u[2] + u[3] - u[7] + u[8] - u[9] + 2 u[11] - 2 u[12] +
2 u[13] + 3 u[14] - 3 u[15] - 2 u[16] + 2 u[17] - 2 u[18] +
3 u[19])} *)

SolveAlways with a arbitrary variable r gives you the desired result



SolveAlways[Thread[coeu == 0], r] // Timing

(* {0.015, {{u[2] -> -u[4] - u[6] - u[8] - 2 u[12] + 2 u[13] - 2 u[15] -
u[16] + u[17] - 3 u[18] + 2 u[19],
u[5] -> u[4] + u[6] - u[10] + u[14] - u[15] + u[19],
u[7] -> -u[10] - u[13],
u[9] -> 2/3 + u[8] + u[10] + u[12] - u[13] + u[15] - u[17] +
2 u[18] - 2 u[19],
u[11] ->
1/3 + u[12] - u[13] - u[14] + u[15] + u[16] - u[17] + u[18] -

u[19], u[1] -> -u[6] - u[8] - u[12] - u[14] - u[18],
u[3] -> -u[4] - u[16] - u[19]}}} *)

Despite expectation Solve doesn't give the minimal solution


Solve[Thread[coeu == 0], Array[u, 19]] // Timing

(* {0., {{u[10] -> -(2/3) - u[1] + u[2] - u[3] + u[4] - u[5] + u[6] +
u[7] - u[8] + u[9],
u[12] -> -(1/3) - 2 u[1] + u[2] + 2 u[4] - u[5] + u[7] - u[8] +
u[11], u[13] ->

2/3 + u[1] - u[2] + u[3] - u[4] + u[5] - u[6] - 2 u[7] + u[8] -
u[9], u[16] ->
2/3 + u[1] - u[2] - u[4] - u[7] + u[8] - u[9] + u[14] - u[15],
u[17] -> 1 - u[2] - u[5] - u[9] - u[11] - u[15],
u[18] ->
1/3 + u[1] - u[2] - 2 u[4] + u[5] - u[6] - u[7] - u[11] - u[14],
u[19] -> -(2/3) - u[1] + u[2] - u[3] + u[7] - u[8] + u[9] - u[14] +
u[15]}}} *)

But Solve with reduced variable set does like SolveAlways.



Solve[Thread[coeu == 0], Array[u, 12]] // Timing

Comments

Popular posts from this blog

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...