# Very slow solve and problem when solving in a nested for loop

I am using a nested for loop to solve a system of equation for {a, ga} for 2 varying parameters: here it's rm and sig.

 params[ss_, mm_] := {c -> 0.05, mu -> 1, bo -> 10, k -> 0.01,  sig -> ss, eps -> 1, rm -> mm}
b[a_] := bo a/(1 + a);
rS[ga_] :=  rm/(1 + c ga);
rI[ga_] := eps rS[ga] ;
getres[aa_, gg_] := Solve[{0 == (rS[gg] x + rI[gg] y) (1 - k (x + y)) - (mu + (b[aa] y)) x + gg y,
0 == b[aa] x y - (mu + aa + gg ) y}, {x, y}][[4]]

par[a_, ga_] :=  Block[{}, res = getres[a, ga]; Seq = x /. res; Ieq = y /. res;
b'[a] - b[a]/(mu + a + ga + sig b[a] Ieq)]

hit[a_, ga_] :=  Block[{}, res = getres[a, ga]; Seq = x /. res; Ieq = y /. res;
rS'[ga]  (1 - k (Seq + Ieq)) + (rI'[ga] (1 - k (Seq + Ieq)) b[a] Ieq)/(mu + a + ga) + (rS[ga] (1 - k (Seq + Ieq)) - mu)/(mu + a + ga)]

getsol[sig_, rm_] :=  Solve[{0 == par[a, ga], 0 == hit[a, ga]} //. params[sig, rm], {a, ga}]


When I write the parameter sig and rm that I want by hand, it is slow but ok:

    getsol[0, 1]

(*{{a -> -0.274696 + 6.20784 I,
ga -> -39.4619 - 3.41054 I}, {a -> -0.274696 - 6.20784 I,
ga -> -39.4619 + 3.41054 I}, {a -> -7.28491*10^-10,
ga -> -1.}, {a -> 7.28491*10^-10, ga -> -1.}, {a -> 0.0358705,
ga -> -0.998713}, {a -> 0.453121, ga -> -0.794681}}*)


but if I created a nested for loop to vary those parameter, it does not solve anything (I know that for some parameter values there is no solution, but for some values, 1 or 2 exist, so I should find something !)

    For [i = 0, i < 0.3, 0.1,
For[j = 0, j < 0.3, 0.1, sig = i; rm = j;
getsol[sig, rm]
]
]


It gives me the error: (Part::partw: Part 4 of {{x->((1+a) (a+ga+mu))/(a bo),y->-(((1+a) mu (a+ga+mu))/(a bo (a+mu)))},{x->0,y->0}} does not exist. >>)

So I know that I'm taking the 4th solution of 'getres' (positive and real), but why doesn't it work here? I've been adding Clear[sig, rm] everywhere but it does not solve the problem ! Any kind of help to improve the time to process and to understand why it does not solve in my for loop is highly welcome !!!

• The first term evaluated is getfol[0,0] and this also gives the same error. So the loop is a red herring. Jan 29, 2015 at 16:00

As noted by Bill S, getsol[0,0] fails. In fact, getsol[sig,0] appears to fail for any value of sig. I suggest you try

For[i = 0, i < 0.3, i = i + 0.1,
For[j = .1, j < 0.3, j = j + 0.1, sig = i; rm = j;
Print[sig, "  ", rm, "  ", getsol[sig, rm]]]]


which excludes rm = 0. Note also that the third argument of the outer For has been changed from 0.1 to i = i + 0.1, and similarly for the inner For. Otherwise, the loops compute the same expression repeatedly. Additionally, Print has been added in order to display the results.

As you remarked, the calculations are very slow for some parameters.

• Awesome !! It works !! Thank you ! Jan 29, 2015 at 17:07

Unnecessary use of SetDelayed rather than Set causes repetitive recalculations and slows everything down. Use of Simplify with some of the interim calculations also improves performance.

params[ss_, mm_] = {c -> 0.05, mu -> 1, bo -> 10, k -> 0.01, sig -> ss,
eps -> 1, rm -> mm};

b[a_] = bo a/(1 + a);

rS[ga_] = rm/(1 + c ga);

rI[ga_] = eps rS[ga];

getres[aa_, gg_] =
Solve[{0 == (rS[gg] x + rI[gg] y) (1 - k (x + y)) - (mu + (b[aa] y)) x +
gg y, 0 == b[aa] x y - (mu + aa + gg) y}, {x, y}][[4]] // Simplify;

par[a_, ga_] = Block[{},
res = getres[a, ga]; Seq = x /. res; Ieq = y /. res;
b'[a] - b[a]/(mu + a + ga + sig b[a] Ieq)];

hit[a_, ga_] = Block[{},
res = getres[a, ga]; Seq = x /. res; Ieq = y /. res;
rS'[ga] (1 - k (Seq + Ieq)) + (rI'[ga] (1 - k (Seq + Ieq)) b[a] Ieq)/(mu +
a + ga) + (rS[ga] (1 - k (Seq + Ieq)) - mu)/(mu + a + ga)] //
Simplify;

getsol[sig_, rm_] :=
Solve[{0 == par[a, ga], 0 == hit[a, ga]} /. params[sig, rm], {a, ga}]

Table[{
"sig = " <> ToString[sig],
"rm = " <> ToString[rm],
getsol[sig, rm]},
{sig, 0, 0.3, 0.1},
{rm, 0.1, 0.3, 0.1}]


{{{"sig = 0.", "rm = 0.1", {{a -> -0.0250166 + 5.01892 I, ga -> -26.1889 - 0.251113 I}, {a -> -0.0250166 - 5.01892 I, ga -> -26.1889 + 0.251113 I}, {a -> -0.901488, ga -> -0.187319}}}, {"sig = 0.", "rm = 0.2", {{a -> -0.0514703 + 5.26649 I, ga -> -28.7332 - 0.542135 I}, {a -> -0.0514703 - 5.26649 I, ga -> -28.7332 + 0.542135 I}, {a -> -0.801944, ga -> -0.356886}}}, {"sig = 0.", "rm = 0.3", {{a -> -0.078452 + 5.4485 I, ga -> -30.68 - 0.854892 I}, {a -> -0.078452 - 5.4485 I, ga -> -30.68 + 0.854892 I}, {a -> -0.701347, ga -> -0.508113}}}}, {{"sig = 0.1", "rm = 0.1", {{a -> -0.0916848 + 5.01913 I, ga -> -26.2134 - 0.256992 I}, {a -> -0.0916848 - 5.01913 I, ga -> -26.2134 + 0.256992 I}, {a -> -0.84974 + 0.0825409 I, ga -> -9.93676 - 15.6267 I}, {a -> -0.84974 - 0.0825409 I, ga -> -9.93676 + 15.6267 I}}}, {"sig = 0.1", "rm = 0.2", {{a -> -0.125137 + 5.26702 I, ga -> -28.7739 - 0.551852 I}, {a -> -0.125137 - 5.26702 I, ga -> -28.7739 + 0.551852 I}, {a -> -0.799511 + 0.0930563 I, ga -> -0.654073 - 16.856 I}, {a -> -0.799511 - 0.0930563 I, ga -> -0.654073 + 16.856 I}}}, {"sig = 0.1", "rm = 0.3", {{a -> -0.157543 + 5.44958 I, ga -> -30.7361 - 0.867882 I}, {a -> -0.157543 - 5.44958 I, ga -> -30.7361 + 0.867882 I}, {a -> -0.748851 + 0.0735576 I, ga -> 7.93714 - 12.7668 I}, {a -> -0.748851 - 0.0735576 I, ga -> 7.93714 + 12.7668 I}}}}, {{"sig = 0.2", "rm = 0.1", {{a -> -0.158886 + 5.02191 I, ga -> -26.2384 - 0.261619 I}, {a -> -0.158886 - 5.02191 I, ga -> -26.2384 + 0.261619 I}}}, {"sig = 0.2", "rm = 0.2", {{a -> -0.199348 + 5.27057 I, ga -> -28.8155 - 0.559327 I}, {a -> -0.199348 - 5.27057 I, ga -> -28.8155 + 0.559327 I}}}, {"sig = 0.2", "rm = 0.3", {{a -> -0.237147 + 5.45404 I, ga -> -30.7934 - 0.877603 I}, {a -> -0.237147 - 5.45404 I, ga -> -30.7934 + 0.877603 I}, {a -> -0.646544 + 0.164275 I, ga -> -5.3496 - 14.9584 I}, {a -> -0.646544 - 0.164275 I, ga -> -5.3496 + 14.9584 I}}}}, {{"sig = 0.3", "rm = 0.1", {{a -> -0.226353 + 5.02736 I, ga -> -26.2639 - 0.264911 I}, {a -> -0.226353 - 5.02736 I, ga -> -26.2639 + 0.264911 I}}}, {"sig = 0.3", "rm = 0.2", {{a -> -0.273768 + 5.27722 I, ga -> -28.8578 - 0.564416 I}, {a -> -0.273768 - 5.27722 I, ga -> -28.8578 + 0.564416 I}}}, {"sig = 0.3", "rm = 0.3", {{a -> -0.316875 + 5.46196 I, ga -> -30.8516 - 0.88386 I}, {a -> -0.316875 - 5.46196 I, ga -> -30.8516 + 0.88386 I}}}}}

• This is a huge improvement to my code and time process ! I always thought SetDelayed was IMPROVING my time, not slowing it down ! Now what took 400s takes 18s ! Thank you ! Jan 30, 2015 at 15:22