See

• The definition of the ODE solver return value ode.h
• The description ODE

<project-file type=“source”/>

<content> #include <gsl/gsl_errno.h> #include <gsl/gsl_matrix.h> #include <gsl/gsl_odeiv.h> #include <math.h> #include <vector> #include <iostream>

#include “ode.h”

ODEResults::ODEResults(int max_household_size, int max_external_degree) {

  MaxHouseholdSize = max_household_size;
  MaxExternalDegree = max_external_degree;

}

void ODEResults::addPoint(double t, double y[]) {

  int mh = MaxHouseholdSize * (MaxHouseholdSize + 3) / 2 + 1;
  int msi = MaxExternalDegree * (MaxExternalDegree + 3) / 2 + 1;
  int m = mh + msi * 2;

  int last = T.size();
  // copy T
  T.push_back(t);
  // copy H
  H.push_back(std::vector<float>());
  H[last].reserve(mh);
  H[last].resize(mh);
  for (int i = 0; i < mh; i++) H[last][i] = y[i];
  //copy S
  S.push_back(std::vector<float>());
  S[last].reserve(msi);
  S[last].resize(msi);
  for (int i = 0; i < msi; i++) S[last][i] = y[i + mh];
  // copy I
  I.push_back(std::vector<float>());
  I[last].reserve(msi);
  I[last].resize(msi);
  for (int i = 0; i < msi; i++) I[last][i] = y[i + mh + msi];

}

float ODEResults::It(int i) const {

  float total = 0;
  int n = I[i].size();
  for (int j = 0; j < n; j++) total += I[i][j];
  return total;

}

float ODEResults::IH(int i) const {

  float total = 0;
  for (int household_size = 1; household_size <= MaxHouseholdSize; household_size++)
      for (int j = 0; j <= household_size; j++)
          total += H[i][ke_index(household_size, j)] * j;
  return total;

}

ODEResults *ode_simulation(ODEParameters &p, double t0, const NetworkStructure &net) {

  double h = 1e-2;
  double epsilon = 1e-2;
  int status;
  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 gsl_odeiv_step_type * T = gsl_odeiv_step_rkf45;
  gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, m);
  gsl_odeiv_control * c = gsl_odeiv_control_y_new (epsilon, 0);
  gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (m);

  gsl_odeiv_system sys = {ode, NULL, m, &p};

  double y[m];

  for (int i = 0; i < mh; i++) y[i] = net.H[i];
  for (int i = 0; i < msi; i++) y[i + mh] = net.S[i];
  for (int i = 0; i < msi; i++) y[i + mh + msi] = net.I[i];

  int lh = net.H.size();
  int ls = net.S.size();
  int li = net.I.size();
  ODEResults *results = new ODEResults(Mh, Msi);
  results->addPoint(t0, y);

  int step = 0;
  while (true) {
      double t = t0 + 1;
      while (t0 < t) {
          status = gsl_odeiv_evolve_apply (e, c, s, &sys, &t0, t, &h, y);
          if (h < 1e-6) h = 1e-6;
          if (status != GSL_SUCCESS) break;
      }
      if (status != GSL_SUCCESS) break;
      ++step;
      results->addPoint(t, y);
      t0 = t;
      if (results->It(step) < 0.5) break;
      if (isnan(results->It(step))) break;
  }

  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (s);
  return results;

} </content>

<use name=“ode.h”/>