Search code examples
c++differential-equations

Euler Method for system of differential equations


I'm trying to recreate computation of a SIR model as described here, with extra midpoint calculations. But for some reason no values actually change during the Euler calculations.

#include <iostream>

using namespace std;

int main(){
    double s[100];
    double i[100];
    double r[100];

    double bb = 1/2;
    double kk = 1/3;
    s[0] = 1;
    i[0] = 1.27e-6;
    r[0] = 0;
    int h = 1;
    int max = 20;

    for (int j = 0; j < max; j += h){
        s[j + 1] = s[j] - h*(bb * s[j] * i[j]);
        i[j + 1] = i[j] + h*(bb * s[j] * i[j] - kk * i[j]);
        r[j + 1] = r[j] + h*(kk * i[j]);

        cout << j << "\t" << s[j] << "\t" << i[j] << "\t" << r[j] << "\n";
        /*
        s[j + 1] = s[j] - 0.5*h*(bb * s[j] * i[j] + bb * s[j+1] * i[j+1]);
        i[j + 1] = i[j] + 0.5*h*(bb * s[j] * i[j] - kk * i[j] + bb * s[j+1] * i[j+1] - kk * i[j+1]);
        r[j + 1] = r[j] + 0.5*h*(kk * i[j] + kk * i[j+1]);
        */
     }
}

Solution

  • Your variables bb and kk are both zero due to integer division. Always use double literals:

    double bb = 1.0/2.0;
    double kk = 1.0/3.0;