1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
| #include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<iostream>
#define dt = 0.000002//integration step
int main();
void EulerStep(double*,double*,double*,double*,double*,double*);
void EulerStep(double *s_, double *i_, double *r_, double *sh_, double *ih_, double *rh_ )
{
/* Ici on doit mettre en place les differentes constantes */
const double pi = 150 ,mu = 0.2,gamma = 0.2, eta = 0.4, omega = 0.6, beta = 0.2, nu = 0.3 ,betar = 0.1 ;
double s, i, r, sh, ih, rh;
double si, ii, ri, shi, rhi, ihi;
double sn, in, rn, shn, rhn, ihn;
s = *s_;i = *i_;r = *r_; sh = *sh_; ih = *ih_; rh = *rh_;
double N = sh+rh+ih;
/* integration at the first order*/
si=s+dt*(pi-(mu+eta)*s+omega*sh);
ii = i + dt*(-(mu+gamma+eta+nu)*i + omega*ih);
ri = r + dt*(nu*i -(mu+eta)*r + omega*rh);
shi = sh + dt*(eta*s -(omega + mu + beta*ih/N)*sh);
ihi = ih + dt*(eta*i +(betar*rh/N + beta*sh/N -(mu + gamma + omega + nu))*ih);
rhi = rh + dt*(nu*ih + eta*r -(omega + mu +betar*ih/N)*rh);
// integration at the second order
sn = s + dt*(pi - (mu + eta)*si + omega*shi);
in = i + dt*(-(mu+gamma+eta+nu)*ii + omega*ihi);
rn = r + dt*(nu*ii -(mu+eta)*ri + omega*rhi);
shn = sh + dt*(eta*si -(omega + mu + beta*ihi/N)*shi);
ihn = ih + dt*(eta*ii +(betar*rhi/N + beta*shi/N -(mu + gamma + omega + nu))*ihi);
rhn = rh + dt*(nu*ihi + eta*ri -(omega + mu +betar*ihi/N)*rhi);
s = (si+sn)/2 ; i = (ii+in)/2 ; r = (ri+rn)/2 ; sh = (shi+shn)/2 ; ih = (ihi+ihn)/2 ; rh = (rhi+rhn)/2 ;
*s_ = s; *i_ = i ; *r_ = r ; *sh_ = sh ; *ih_ = ih ; *rh_ = rh ;
}; |
Partager