0
$\begingroup$

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 !!!

$\endgroup$
  • 2
    $\begingroup$ The first term evaluated is getfol[0,0] and this also gives the same error. So the loop is a red herring. $\endgroup$ – bill s Jan 29 '15 at 16:00
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Awesome !! It works !! Thank you ! $\endgroup$ – user3767071 Jan 29 '15 at 17:07
1
$\begingroup$

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}}}}}

$\endgroup$
  • $\begingroup$ 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 ! $\endgroup$ – user3767071 Jan 30 '15 at 15:22

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.