The code below models a simple SIR model (used in disease control) in Mathematica. (I copied it directly from my notebook).

The equations can be solved using NDSolve and the solutions are inserted into three different functions for further use.

As can be seen the Beta term on the first line varies depending on the value of Inf[t], which is one of the three solutions of the NDSolve function.

This code works fine and I have included this in order to better explain my quesion below.

Beta = Piecewise[{{0.01, Inf[t] > 20}, {.06, Inf[t] <= 20}}];
Mu = 0.1;
Pop = 100;
ans = NDSolve[{S'[t] == -Beta S[t] Inf[t], 
    Inf'[t] == Beta S[t] Inf[t] - Mu Inf[t], 
    R'[t] == Mu Inf[t], 
    S[0] == Pop - 1, Inf[0] == 1, 
    R[0] == 0}, {S[t], Inf[t], R[t]}, {t, 0, 10}];
Sus[t_] = S[t] /. ans[[1, 1]];
Infected[t_] = Inf[t] /. ans[[1, 2]];
Rec[t_] = R[t] /. ans[[1, 3]];

I now wanted to update the code so that instead of having an either/or value for the Beta parameter based on the Inf[t] value, I would have the Beta value being equal to the output of a function where Inf[t] is the input. This can be seen below where UpdateTransmission[] is the function.

When I try and run the code below though the Beta value remains at 0 and does not increase. The problem is not with the UpdateTransmission function as I have tested this independently.

Beta = UpdateTransmission[SpinMatrix, ThresholdMatrix, Inf[t]];
Mu = 0.1;
Pop = 100;
ans = NDSolve[{S'[t] == -Beta S[t] Inf[t], 
    Inf'[t] == Beta S[t] Inf[t] - Mu I开发者_StackOverflow中文版nf[t], 
    R'[t] == Mu Inf[t], S[0] == Pop - 1, Inf[0] == 1, 
    R[0] == 0}, 
    {S[t], Inf[t], R[t]}, {t, 0, 10}];
Sus[t_] = S[t] /. ans[[1, 1]];
Infected[t_] = Inf[t] /. ans[[1, 2]];
Rec[t_] = R[t] /. ans[[1, 3]];

Plot[{Sus[t], Infected[t], Rec[t]}, {t, 0, 5}]

Can anyone shed some light on why this may not be running correctly?

Edit: here is the updated function

UpdateTransmission[S_, Th_, Infect_] := Module[{BetaOverall},
 P = S;
 For[i = 1, i <= Pop, i++,
    P[[i]] = Sign[Infect - Th[[i]]];];
   BetaOverall = ((Count[P, 1]*.02) + (Count[P, -1]*.5))/Pop
]

Here are the two lists that are referred to in the code above:

SpinMatrix = Table[-1, {Pop}]

val := RandomReal[NormalDistribution[.5, .1]]
ThresholdMatrix = Table[Pop*val, {Pop}]

Edit Edit

Ok I've put everything together and tried to plot my three curves. Now as can be seen here they are all flat-lining. The Sus[t] line is staying at 100 whilst the other two seem to be staying below 1. What should be happening here is that the Sus[t] line should drop considerably and the other two lines should ramp up.

(I tried to insert and image but I can't as I don't have the reputation points required so I'll just past in the code and you can see the plot yourself on your own machine)

 Pop = 100;
SpinMatrix = Table[-1, {Pop}];
val := RandomReal[NormalDistribution[.5, .1]];
ThresholdMatrix = Table[Pop*val, {Pop}];

updateTransmission[S_, Th_, Infect_] := Module[{}, P = S;
   For[i = 1, i <= Pop, i++, P[[i]] = Sign[Infect - Th[[i]]];];
   Return[((Count[P, 1]*.02) + (Count[P, -1]*.5))/Pop]];

beta[t_] := updateTransmission[SpinMatrix, ThresholdMatrix, Inf[t]];
mu = 0.1;
ans = NDSolve[{S'[t] == -beta[t] S[t] Inf[t], 
    Inf'[t] == beta[t] S[t] Inf[t] -
      mu Inf[t], R'[t] == mu Inf[t], S[0] == Pop - 1, Inf[0] == 1, 
    R[0] == 0}, {S[t], Inf[t], R[t]}, {t, 0, 10}];
Sus[t_] = S[t] /. First@ans;
Infected[t_] = Inf[t] /. First@ans;
Rec[t_] = R[t] /. First@ans;
Plot[{Sus[t], Infected[t], Rec[t]}, {t, 0, 10}]

The output that I am expecting should look similar to that of the code given below:

Beta = Piecewise[{{0.5, Inf[t] > 20}, {.02, Inf[t] <= 20}}];
Mu = 0.1;
Pop = 100;
ans = NDSolve[{S'[t] == -Beta S[t] Inf[t], 
    Inf'[t] == Beta S[t] Inf[t] - Mu Inf[t], 
    R'[t] == Mu Inf[t], S[0] == Pop - 1, Inf[0] == 1, 
    R[0] == 0}, {S[t], Inf[t], R[t]}, {t, 0, 10}];
Sus[t_] = S[t] /. ans[[1, 1]];
Infected[t_] = Inf[t] /. ans[[1, 2]];
Rec[t_] = R[t] /. ans[[1, 3]];
Plot[{Sus[t], Infected[t], Rec[t]}, {t, 0, 10}]

The culprit is Sign[ ]

I don't know why, but I traced the problem to the Sign[ ] function that is not working properly inside NDSolve!

Removing it:

Pop = 100;
SpinMatrix = Table[-1, {Pop}];
val := RandomReal[NormalDistribution[.5, .1]];

ThresholdMatrix = Table[Pop*val, {Pop}];

updateTransmission[Th_, Inf_] :=
  Total[Table[If[Inf >= Th[[i]], 2/100, 1/2]/Pop , {i, Pop}]];

beta[t_] := updateTransmission[ThresholdMatrix, Inf[t]];
mu = 0.1;

ans = NDSolve[{
    S'[t] == -beta[t] S[t] Inf[t],
    Inf'[t] == beta[t] S[t] Inf[t] - mu Inf[t],
    R'[t] == mu Inf[t],
    S[0] == Pop - 1,
    R[0] == 0,
    Inf[0] == 1}, {S[t], Inf[t], R[t]}, {t, 0, 10}];

Sus[t_] = S[t] /. First@ans;
Infected[t_] = Inf[t] /. First@ans;
Rec[t_] := R[t] /. First@ans;
Plot[{Sus[t], Infected[t], Rec[t]}, {t, 0, 10}]  

Gives:

Probably someone with better knowledge of Mma could explain what is happening in your code.

HTH!

In some ways, you are encountering the difference between Set (=) and SetDelayed (:=). For instance, if you wrote f = 7, f becomes 7 in all occurrences of f after it was initialized. But, if you wrote f = 7 t instead, and tried to use it as you would a function, i.e. f[8], you'd get (7 t)[8] because Set says that the value of f is unchanging. SetDelayed, however, implies that the value of f will change and must be reevaluated every time it occurs. Your initial case, though, is special.

When you wrote

Beta = Piecewise[{{0.01, Inf[t] > 20}, {.06, Inf[t] <= 20}}]

Inf[t] was undefined, so that it remained unevaluated. But, every occurence of Beta in your differential equations was replaced by the above formula, courtesy of Set, so NDSolve only saw the Piecewise functions. In your second case, you wrote

Beta = UpdateTransmission[Inf[t]]

Here the problem is that UpdateTransmission is executed only when Beta is initially defined, and while Piecewise remains unevaluated, UpdateTransmission most likely still gives a result for a purely symbolic input. I'd try one of three things,

1. replace every occurrence of Beta in you equations with UpdateTransmission[Inf[t]],

2. redefine Beta using SetDelayed, e.g.

   Beta := UpdateTransmission[Inf[t]]
   

   so that it will be reevaluated every time it is encountered, or

3. redefine UpdateTransmission to not accept symbols via either

   UpdateTransmission[x_?(Head[#]=!=Symbol&)] := ...
   

   or

   UpdateTransmission[x_] /; Head[x]=!= Symbol := ...
   

Option 3 works by forcing UpdateTransmission[Inf[t]] to remain unevaluated, and effectively does the same thing as option 1. But, it requires a minimum of change. Personally, I'm in favor of options 1 or 3, as I don't know how many times Beta will need to be reevaluated as NDSolve operates.