Приведенный ниже код моделирует простую модель SIR (используется для контроля заболеваний) в Mathematica. (Я скопировал его прямо из блокнота).
Уравнения могут быть решены с помощью NDSolve, а решения вставлены в три разные функции для дальнейшего использования.
Как видно, термин Бета в первой строке изменяется в зависимости от значения Inf [t], которое является одним из трех решений функции NDSolve.
Этот код работает нормально, и я включил его, чтобы лучше объяснить мой вопрос ниже.
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]];
Теперь я хотел обновить код таким образом, чтобы вместо наличия либо / или значения для параметра Beta на основе значения Inf[t], у меня было бы значение Beta, равное выходу функции где Inf[t] является входом. Это можно увидеть ниже, где UpdateTransmission[] - функция.
Когда я пытаюсь запустить приведенный ниже код, значение Beta остается равным 0 и не увеличивается. Проблема не в функции UpdateTransmission, поскольку я проверил это независимо.
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 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, 5}]
Может кто-нибудь пролить свет на то, почему это может работать неправильно?
Редактировать: вот обновленная функция
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
]
Вот два списка, которые упоминаются в приведенном выше коде:
SpinMatrix = Table[-1, {Pop}]
val := RandomReal[NormalDistribution[.5, .1]]
ThresholdMatrix = Table[Pop*val, {Pop}]
Редактировать Редактировать
Хорошо, я собрал все вместе и попытался построить три кривые. Теперь, как видно здесь, все они плоские. Линия Sus [t] остается на уровне 100, в то время как две другие, похоже, остаются ниже 1. Что здесь должно произойти, так это то, что линия Sus [t] должна значительно упасть, а две другие линии должны увеличиться.
(Я пытался вставить и изобразить, но не могу, поскольку у меня нет требуемых очков репутации, поэтому я просто пройду мимо кода и вы сами сможете увидеть сюжет на своей машине)
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}]
Вывод, который я ожидаю, должен выглядеть примерно так, как показано в приведенном ниже коде:
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}]