See
- The definition of the ODE parameters equations.h
- The description Equations
<project-file type=“source”/> <content> #include “equations.h” #include <gsl/gsl_errno.h> #include “ke_index.h”
This function defines the index of the variable when its degree is k, and the infectious degree (we used to call it the effective degree) is e. /* inline int ke_index(int k,int e) { return 1)
1)
k-1)*(k-1+3)/2 + e + 1);
}
*/
inline double V(const double X[], int n, int i, int max)
{
if (n < 0) return 0; if (n > max) return 0; if (i > n) return 0; if (i < 0) return 0; return X[ke_index(n, i)];} int ode(double t, const double y[], double dydt[], void *params) {
const ODEParameters *p = (ODEParameters*) params; int Mh = p->MaxHouseholdSize; int mh = Mh * (Mh + 3) / 2 + 1; int Msi = p->MaxExternalDegree; int msi = Msi * (Msi + 3) / 2 + 1; int m = mh + msi * 2;
const double *H = y, *S = &y[mh], *I = &y[mh + msi]; double beta_h = p->BetaHousehold; double beta_e = p->BetaExternal; double gamma = p->Gamma;
double external_infection_rate = 0; for (int k = 0; k <= Msi; k++) for (int i = 0; i <= k; i++) external_infection_rate += beta_e * i * V(S, k, i, Msi); double household_susceptible = 0; for (int n = 0; n <= Mh; n++) for (int j = 0; j <= n; j++) household_susceptible += pnj(n, j) * V(H, n, j, Mh); // the H equations for (int n = 0; n <= Mh; n++) for (int j = 0; j <= n; j++) { dydt[ke_index(n, j)] = beta_h * ((n - j + 1) * (j - 1) * V(H, n, j - 1, Mh) - (n - j) * j * V(H, n, j, Mh)) + gamma * ((j + 1) * V(H, n + 1, j + 1, Mh) - j * V(H, n, j, Mh)) + external_infection_rate / household_susceptible * (pnj(n, j - 1) * V(H, n, j - 1, Mh) - pnj(n, j) * V(H, n, j, Mh)); } double incidence_s_neighbors = 0; for (int k = 0; k <= Msi; k++) for (int i = 0; i <= k; i++) { incidence_s_neighbors += beta_e * (k - i) * i * V(S, k, i, Msi); } double susceptible_s_neighbors = 0; for (int k = 0; k <= Msi; k++) for (int i = 0; i <= k; i++) { susceptible_s_neighbors += (k - i) * V(S, k, i, Msi); } double internal_infection_rate = 0; for (int n = 0; n <= Mh; n++) for (int j = 0; j <= n; j++) internal_infection_rate += beta_h * (n - j) * j * V(H, n, j, Mh); double all_s = 0; for (int j = 0; j < msi; j++) all_s += S[j]; // The S equations for (int k = 0; k <= Msi; k++) for (int i = 0; i <= k; i++) { int s = k - i; dydt[mh + ke_index(k, i)] = - beta_e * i * V(S, k, i, Msi) + gamma * ((i + 1) * V(S, k + 1, i + 1, Msi) - i * V(S, k, i, Msi)) + incidence_s_neighbors / susceptible_s_neighbors * ((s + 1) * V(S, k, i - 1, Msi) - s * V(S, k, i, Msi)) - internal_infection_rate / all_s * V(S, k, i, Msi) + internal_infection_rate / all_s * ((s + 1) * V(S, k, i - 1, Msi) - s * V(S, k, i, Msi)); } double incidence_i_neighbors = 0; for (int k = 0; k <= Msi; k++) for (int i = 0; i <= k; i++) { incidence_i_neighbors += beta_e * i * i * V(S, k, i, Msi); } double infectious_s_neighbors = 0; for (int k = 0; k <= Msi; k++) for (int i = 0; i <= k; i++) { infectious_s_neighbors += (k - i) * V(I, k, i, Msi); } // The I equations for (int k = 0; k <= Msi; k++) for (int i = 0; i <= k; i++) { int s = k - i; dydt[mh + msi + ke_index(k, i)] = beta_e * i * V(S, k, i, Msi) - gamma * V(I, k, i, Msi) + gamma * ((i + 1) * V(I, k + 1, i + 1, Msi) - i * V(I, k, i, Msi)) + incidence_i_neighbors / infectious_s_neighbors * ((s + 1) * V(I, k, i - 1, Msi) - s * V(I, k, i, Msi)) + internal_infection_rate / all_s * V(S, k, i, Msi) + internal_infection_rate / all_s * ((s + 1) * V(I, k, i - 1, Msi) - s * V(I, k, i, Msi)); } return GSL_SUCCESS;} </content> <use name=“equations.h”/> <use name=“ke_index.h”/>