# speeding up Solve or NSolve

I have a system of equations which I am trying to solve.

Solve[ms^2 == -2 (M^2 - (6 (k1 + k2) b^2 m^2)/(c + M^2 + 2 m5^2)^2 +
c + 2 m5^2) &&
ma^2 == -2 (M^2 - (2 (3 k1 + k2) b^2 m^2)/(c + M^2 + 2 m5^2)^2 -
c + 2 m5^2) &&
mp^2 == 2 b m (Sqrt[(c + M^2 + 2 M5^2)/(2 (k1 + k2))] +
m b/(2 (c + M^2 + 2 M5^2)))^(-1) &&
v == Sqrt[(c + M^2 + 2 M5^2)/(2 (k1 + k2))] +
m b/(2 (c + M^2 + 2 M5^2)), {k1, k2, c, b}] /. {ms -> 500,
ma -> 980, v -> 92, m -> 5.5, M -> 364, m5 -> 0}


But wolfram was trying to solve it more than 3 hours. Could you tell me what I'm doing wrong?

• If you want to get anything useful from NSolve it will require that there be as many equations as unknowns. Right now it is too few (4 vs 6). Also the substitutions should happen before calling Solve, in order for there to be much chance it will run to completion. Commented Mar 8, 2017 at 21:33
• Do you know M5 and m5 are different variables.? I notice you don't give value to M5 when you are done. Commented Mar 9, 2017 at 16:41
• Yes, I know about it, but let M5 and m5 are same variables. Commented Mar 10, 2017 at 17:23

One approach would be to convert to algebraic expressions, do the numeric substitutions, then find a Groebner basis that in effect "triangularizes" it so that the variables of interest can be solved in terms of the remaining parameters.

For this purpose I show a direct code below. One could explicitly change the radicals into new variables, with corresponding relations, but I do not go that far. GroebnerBasis will figure it out anyway, and the only drawback is once it removes its "internal" variables from the output, things involving the radicals might look weird.

  eqns = Apply[List,
ms^2 == -2 (M^2 - (6 (k1 + k2) b^2 m^2)/(c + M^2 + 2 m5^2)^2 +
c + 2 m5^2) &&
ma^2 == -2 (M^2 - (2 (3 k1 + k2) b^2 m^2)/(c + M^2 + 2 m5^2)^2 -
c + 2 m5^2) &&
mp^2 == 2 b m (Sqrt[(c + M^2 + 2 M5^2)/(2 (k1 + k2))] +
m b/(2 (c + M^2 + 2 M5^2)))^(-1) &&
v == Sqrt[(c + M^2 + 2 M5^2)/(2 (k1 + k2))] +
m b/(2 (c + M^2 + 2 M5^2))] /. {ms -> 500, ma -> 980, v -> 92,
m -> 11/2, M -> 364, m5 -> 0};
exprs = Apply[Subtract, eqns, {1}]
vars = Variables[exprs]

Timing[
gbrat = GroebnerBasis[exprs, {k1, k2, c, b},
CoefficientDomain -> RationalFunctions];]

(* Out[71]= {0.742745, Null} *)


The result is on the long side.

In[72]:= gbrat

(* Out[72]= {-11 b + 92 mp^2,
1269701225750632387225583616 + 16 c^5 +
38331760226742917136384 M5^2 + 289305037335035904 M5^4 -
4791470028342864642048 mp^2 - 72326259333758976 M5^2 mp^2 -
23391518267959296 mp^4 - 1263973681152 M5^2 mp^4 -
19079424 M5^4 mp^4 - 96 M5^6 mp^4 +
c^4 (12599680 + 64 M5^2 - 8 mp^2) +
c^3 (3868798402560 + 41918976 M5^2 + 64 M5^4 - 5239872 mp^2 -
16 M5^2 mp^2 - 11 mp^4) +
c^2 (582821073209589760 + 9921096966144 M5^2 + 33439232 M5^4 -
1240137120768 mp^2 - 8359808 M5^2 mp^2 - 4247368 mp^4 -
72 M5^2 mp^4) +
c (43262715370708774420480 + 1016778629212143616 M5^2 +
5490532483072 M5^4 - 127097328651517952 mp^2 -
1372633120768 M5^2 mp^2 - 546197270528 mp^4 -
19079424 M5^2 mp^4 - 144 M5^4 mp^4),
9582940056685729284096 + 16 c^4 + 144652518667517952 M5^2 -
36163129666879488 mp^2 - 176545090176 mp^4 - 6874800 M5^2 mp^4 -
44 M5^4 mp^4 + c^3 (10479744 + 32 M5^2 - 8 mp^2) +
c^2 (2480274241536 + 16719616 M5^2 - 4179904 mp^2 - 11 mp^4) +
c (254194657303035904 + 2745266241536 M5^2 - 686316560384 mp^2 -
2789912 mp^4 - 50 M5^2 mp^4) + (500000 M5^4 mp^4 - 8 M5^6 mp^4)/(
132496 + c + 2 M5^2),
9582940056685729284096 + 16 c^4 + 289305037335035904 M5^2 +
2183500161024 M5^4 - 36163129666879488 mp^2 -
545875040256 M5^2 mp^2 - 176545090176 mp^4 - 9539712 M5^2 mp^4 -
144 M5^4 mp^4 - (96 M5^6 mp^4)/(132496 + c) +
c^3 (10479744 + 64 M5^2 - 8 mp^2) +
c^2 (2480274241536 + 33439232 M5^2 + 64 M5^4 - 4179904 mp^2 -
16 M5^2 mp^2 - 11 mp^4) +
c (254194657303035904 + 5490532483072 M5^2 + 24959488 M5^4 -
686316560384 mp^2 - 6239872 M5^2 mp^2 - 2789912 mp^4 -
72 M5^2 mp^4),
3117801746841600 + 29507389184 c - 87392 c^2 - c^3 +
4232 k2 mp^4, -6936898224442368 - 87155868800 c - 130156 c^2 +
c^3 + 6348 k1 mp^4, -29703370839766117917696000000000000 +
520066650864685777821696000000 M5^2 +
13622605280123027498139648 M5^4 - 238099500614170484736 M5^6 -
1931758864590848 M5^8 + 17181784576 M5^10 - 67712 M5^12 +
c^4 (-49593750000000 + 2116000000 M5^2 -
25392 M5^4) + (-1953125000000000 M5^8 + 93750000000 M5^10 -
1500000 M5^12 + 8 M5^14)/(k1 + k2) +
112091575744800288000000000000 mp^2 -
3654574303540812288000000 M5^2 mp^2 +
3757441485177747456 M5^4 mp^2 + 841797996443136 M5^6 mp^2 -
5416892288 M5^8 mp^2 + 16928 M5^10 mp^2 +
547220816619750000000000 mp^4 + 3467897167974000000 M5^2 mp^4 -
551724374390688 M5^4 mp^4 + 1122415272 M5^6 mp^4 +
125902 M5^8 mp^4 +
c^3 (-32483112750000000000 + 1220633644000000 M5^2 -
9754353728 M5^4 - 76176 M5^6 + 24796875000000 mp^2 -
1058000000 M5^2 mp^2 + 12696 M5^4 mp^2) +
c^2 (-7687881291636000000000000 + 241642470943136000000 M5^2 -
409170245317632 M5^4 - 37685045888 M5^6 - 16928 M5^8 +
12956069625000000000 mp^2 - 519729804000000 M5^2 mp^2 +
5311007648 M5^4 mp^2 + 12696 M5^6 mp^2 + 34095703125000 mp^4 -
1454750000 M5^2 mp^4 + 17457 M5^4 mp^4) +
c (-787904142851402304000000000000 +
19435155129766498304000000 M5^2 + 160779731080178020352 M5^4 -
5577633215976448 M5^6 - 10833784576 M5^8 + 33856 M5^10 +
2127313244784000000000000 mp^2 -
77871254610784000000 M5^2 mp^2 + 540357461329408 M5^4 mp^2 +
6538338432 M5^6 mp^2 - 16928 M5^8 mp^2 +
8647637390625000000 mp^4 - 168524455750000 M5^2 mp^4 -
4003347156 M5^4 mp^4 + 96807 M5^6 mp^4), -1 +
92/((11 b)/(4 (132496 + c + 2 M5^2)) + Sqrt[(132496 + c + 2 M5^2)/(
k1 + k2)]/Sqrt[2]),
220407621303771773534208 + 368 c^4 + 3327007929352912896 M5^2 -
831751982338228224 mp^2 + 46000000 M5^4 mp^2 - 736 M5^6 mp^2 -
4060537074048 mp^4 - 158120400 M5^2 mp^4 - 1012 M5^4 mp^4 +
c^3 (241034112 + 736 M5^2 - 184 mp^2) +
Sqrt[(132496 + c + 2 M5^2)/(
k1 + k2)] (-250000 Sqrt[2] M5^4 mp^2 + 4 Sqrt[2] M5^6 mp^2) +
c^2 (57046307555328 + 384551168 M5^2 - 96137792 mp^2 - 253 mp^4) +
c (5846477117969825792 + 63141123555328 M5^2 -
15785280888832 mp^2 - 64167976 mp^4 - 1150 M5^2 mp^4)} *)


The point is one can now "back substitute" to get solutions. The first polynomial implies b = 92/11*mp^2, the second one gives c as a quintic in the parameters, etc.

• It is very interesting answer! Thank you! It would be interesting to know is it posible to do for more complicated problem? For example if I try to find a numerical calculation of a integral equation with parameters without FindRoot. Commented Mar 10, 2017 at 17:16
• For example, I know how to solve: Veffd2i[M_?NumberQ, m_?NumberQ, T_?NumberQ] := 4 + NIntegrate[-((2 E^(-((m + Sqrt[4 + (l + M)^2])/T)) l^2)/((1 + E^(-((m + Sqrt[(l + M)^2 + 4])/T))) T Sqrt[ 4 + (l + M)^2])), {l, -100, 100}, Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0}]; if I use FindRoot: Veffd2i2[m_?NumberQ, T_?NumberQ] := M /. FindRoot[Veffd2i[M, m, T] == 0, {M, 0}]; Veffd2i2[1, -2] Commented Mar 10, 2017 at 17:20
• I see no way to extending a method using polynomial algebra to handle root-finding of what amount to "black-box" functions that use e.g. quadrature. Methods using the likes of FindRoot or perhaps path-tracking via NDSolve, are the ways I know for how to go about this. Commented Mar 10, 2017 at 21:07

To be able to get numerical solution. you will need to specify all the parameters. Choosing random values for M5 and mp,

ms = 500; ma = 980; v = 92; m = 5.5; M = 364; m5 = 0; M5 = 10; mp = 10;

Eq1 = ms^2 == -2 (M^2 - (6 (k1 + k2) b^2 m^2)/(c + M^2 + 2 m5^2)^2 + c + 2 m5^2);
Eq2 = ma^2 == -2 (M^2 - (2 (3 k1 + k2) b^2 m^2)/(c + M^2 + 2 m5^2)^2 - c + 2 m5^2);
Eq3 = mp^2 == 2 b m (Sqrt[(c + M^2 + 2 M5^2)/(2 (k1 + k2))] +
m b/(2 (c + M^2 + 2 M5^2)))^(-1);
Eq4 = v == Sqrt[(c + M^2 + 2 M5^2)/(2 (k1 + k2))] + m b/(2 (c + M^2 + 2 M5^2));

NSolve[{Eq1, Eq2, Eq3, Eq4}, {k1, k2, c, b}]


{{k1 -> 1.60642*10^8, k2 -> -1.60642*10^8, c -> -257496., b -> 836.364}, {k1 -> 254.123, k2 -> -224.074, c -> -132671., b -> 836.364}, {k1 -> 0.13213, k2 -> -0.116479, c -> -132492., b -> 836.364}, {k1 -> 0.128457, k2 -> -0.113241, c -> -132500., b -> 836.364}}