gama-platform / gama.old

Main repository for developing the 1.x versions of GAMA
GNU General Public License v3.0
302 stars 99 forks source link

ODE : influence of the integration step #2356

Closed benoitgaudou closed 6 years ago

benoitgaudou commented 6 years ago

Steps to reproduce

  1. Run the model showing the influence of the integration step (e.g. Plugin Models / Ordinary Differential Equations / SIR (Influence of Integration Step) )

Expected behavior

A difference between the 3 integration steps are expected in the result.

Actual behavior

The curves seem identical. If there is no difference, the model seems to be useless.

System and version

MACOSX, GAMA Git Mars, Gama 1.8 preRelease

hqnghi88 commented 6 years ago

I tested wit 1.8 , it seems not identical at all. It faster with smaller h. Did you see this phenomenon? image

benoitgaudou commented 6 years ago

Hi Nghi, You are right on this model, I get the same results. But my point is that that is not the h that has an influence on the results, but the myUnit variable, that changes the equation. Notice that in the model: Library models / Toy Models / Predator Prey / Lotka-Volterra (Influence of Integration Step) the various values of the steps do not have any influence on the result. benoit

hqnghi88 commented 6 years ago

Dear Benoit, Actually it was. There are two thing, unit step (myUnit) of each agent, and the integration step (h) . For myUnit you can see the differ of evolution speed of these charts (is it your remark?). Sorry for the h, there is a hidden information, you can see in the screenshot , i add some "bis" display to show the different of integration step. the smaller step give the finer result (more number values). I did not take in charge the lotka. I review it now. image

hqnghi88 commented 6 years ago

@@ the model to test, i am working on oxygen , no way to safe commit. /**

model SIR_influence_of_integration_step

global { float step<-1#s; float beta <- 0.8 ;
float delta <- 0.01 ;

float s1 <- 1#s;
float s2 <- 1.5#s;
float s3 <- 2#s;

init {
    create SIR_agt with: [h::1,myUnit::s1];
    create SIR_agt with: [h::0.1,myUnit::s2];
    create SIR_agt with: [h::0.01,myUnit::s3];  
}  

}

species SIR_agt { int N <- 500; float t;

float I <- 1.0; 
float S <- N - I; 
float R <- 0.0; 

float h;        
float myUnit<-1#s;
equation SIR{ 
    diff(S,t) =myUnit* (- beta  * S * I / N);
    diff(I,t) =myUnit* (beta * S * I / N) - (delta * I);
    diff(R,t) =myUnit* (delta * I);
} 

reflex solving {
    solve SIR method: "rk4" step: h;
}      

}

experiment mysimulation1 type: gui { output { display SIR_1 { chart "SI - h=1" type: series background: #white {
data 'S' value: first(SIR_agt where (each.myUnit = s1)).S color: #green;
data 'I' value: first(SIR_agt where (each.myUnit = s1)).I color: #red ; data 'R' value: first(SIR_agt where (each.myUnit = s1)).R color: #blue ;
} } display "SIR_1bis" { chart "SI - h=1" type: series background: #white x_serie: first(SIR_agt where (each.myUnit = s1)).t[] {
data 'S' value: first(SIR_agt where (each.myUnit = s1)).S[] color: #green;
data 'I' value: first(SIR_agt where (each.myUnit = s1)).I[] color: #red ; data 'R' value: first(SIR_agt where (each.myUnit = s1)).R[] color: #blue ;
} }

    display "SIR_10" {
        chart "SI - h=0.1" type: series background: #white
        {       
            data 'S' value: first(SIR_agt where (each.myUnit = s2)).S color: #green;                
            data 'I' value: first(SIR_agt where (each.myUnit = s2)).I color: #red ;
            data 'R' value: first(SIR_agt where (each.myUnit = s2)).R color: #blue ;                
        }
    }

    display "SIR_10bis" {
        chart "SI - h=0.1" type: series background: #white
        x_serie:  first(SIR_agt where (each.myUnit = s2)).t[]
        {       
            data 'S' value: first(SIR_agt where (each.myUnit = s2)).S[] color: #green;              
            data 'I' value: first(SIR_agt where (each.myUnit = s2)).I[] color: #red ;
            data 'R' value: first(SIR_agt where (each.myUnit = s2)).R[] color: #blue ;              
        }
    }

    display "SIR_100"  {
        chart "SI - h=0.01" type: series background: #white
        {       
            data 'S' value: first(SIR_agt where (each.myUnit = s3)).S color: #green;                
            data 'I' value: first(SIR_agt where (each.myUnit = s3)).I color: #red ;
            data 'R' value: first(SIR_agt where (each.myUnit = s3)).R color: #blue ;                
        }
    }

    display "SIR_100bis"  {
        chart "SI - h=0.01" type: series background: #white
        x_serie:  first(SIR_agt where (each.myUnit = s3)).t[]
        {       
            data 'S' value: first(SIR_agt where (each.myUnit = s3)).S[] color: #green;              
            data 'I' value: first(SIR_agt where (each.myUnit = s3)).I[] color: #red ;
            data 'R' value: first(SIR_agt where (each.myUnit = s3)).R[] color: #blue ;              
        }
    }       
}

}

hqnghi88 commented 6 years ago

Lotka will have same react if i add myUnit as followed: /**

model ODE_LotkaVolterra_InfluenceTimeStep

global {

float prey_birth_rate ;         // natural birth rate of preys
float predation_rate ;          // death rate of preys due to predators
float predator_death_rate ;     // natural death rate of predators
float predation_efficiency ;    // birth rate of predators due to prey consumption

float nb_prey_init ;            // initial number of preys
float nb_predator_init  ;       // initial number of predators

float integration_step ;    // integration time step of the Runge Kutta 4 method
float t;                    // simulation time : t = n * integration_time_step  where n is the number of already computed time step

float s1 <- 1#s;
float s2 <- 1.5#s;
float s3 <- 2#s;

float integration_time_step1  <- 1.0;  // first integration time step to compare 
float integration_time_step2  <- 0.1;  // second integration time step to compare 
float integration_time_step3  <- 0.01;  // third integration time step to compare 

list<LotkaVolterra_agent> LV_agents;

init{
    create LotkaVolterra_agent number: 1 with:[integration_time_step::integration_time_step1,myUnit::s1];   // creation of an agent containing the ODE model with an integration time step of value integration_time_step1
    create LotkaVolterra_agent number: 1 with:[integration_time_step::integration_time_step2,myUnit::s2];   // creation of an agent containing the ODE model with an integration time step of value integration_time_step2
    create LotkaVolterra_agent number: 1 with:[integration_time_step::integration_time_step3,myUnit::s3];   // creation of an agent containing the ODE model with an integration time step of value integration_time_step3
    LV_agents <- list(LotkaVolterra_agent);
}

}

species LotkaVolterra_agent {

float nb_prey <- nb_prey_init ;                 // number of preys initialized with the values given by the user
float nb_predator <- nb_predator_init ;         // number of predators initialized with the values given by the user

float myUnit<-1#s;

float integration_time_step ;                   // integration time step used in the Runge Kutta 4 method

equation lotka_volterra { 
    diff(nb_prey,t) = myUnit* nb_prey * (prey_birth_rate - predation_rate * nb_predator);                   // evolution of the number of preys duting an integration time step
    diff(nb_predator,t) = myUnit* - nb_predator * (predator_death_rate - predation_efficiency * nb_prey);       // evolution of the number of predator during an integration time step
  }
  reflex solving {        
    solve lotka_volterra method: "rk4" step:integration_time_step;// cycle_length:1/integration_time_step;          // use of runge kutta 4 method with an integration time step of value integration_time_step
   }

}

experiment maths type: gui {

parameter "Prey birth rate" var: prey_birth_rate <- 0.05 min: 0.0 max: 1.0 category: "Prey";                        // the user defines the value of parameter prey_birth_rate on the interface, the default value is 0.05 and this value must be between 0 and 1
parameter "Predation rate" var: predation_rate <- 0.001 min: 0.0 max: 1.0 category: "Prey";                         // the user defines the value of parameter prey_birth_rate on the interface, the default value is 0.001 and this value must be between 0 and 1
parameter "Predator death rate" var: predator_death_rate <- 0.03 min: 0.0 max: 1.0 category: "Predator";            // the user defines the value of parameter predator_death_rate on the interface, the default value is 0.03 and this value must be between 0 and 1
parameter "Predation efficiency" var: predation_efficiency <- 0.0002 min: 0.0 max: 1.0 category: "Predator";        // the user defines the value of parameter predation_efficiency on the interface, the default value is 0.0002 and this value must be between 0 and 1

parameter "Initial number of prey" var: nb_prey_init <- 250.0 min: 1.0 category: "Prey";                            // the user defines the value of parameter predation_efficiency on the interface, the default value is 250, the minimum possible value is 1
parameter "Initial number of predator" var: nb_predator_init <- 45.0 min: 1.0 category: "Predator";                 // the user defines the value of parameter predation_efficiency on the interface, the default value is 45, the minimum possible value is 1

parameter "Integration time step of the first chart " var:  integration_time_step1 <- 1.0  min: 0.0 max:1.0 category: "Integration time steps";     // the user defines the value of the first integration step he wants to compare, the default value is 1 and this value must be between 0 and 1
parameter "Integration time step of the second chart " var:  integration_time_step2 <- 0.1  min: 0.0 max: 1.0 category: "Integration time steps";   // the user defines the value of the second integration step he wants to compare, the default value is 0.1 and this value must be between 0 and 1
parameter "Integration time step of the third chart " var:  integration_time_step3 <- 0.01  min: 0.0 max: 1.0 category: "Integration time steps";   // the user defines the value of the third integration step he wants to compare, the default value is 0.01 and this value must be between 0 and 1

output {        
    display TimeSeries  {   // creation of a display to show time series of the model, values are plotted at every step. Since there is more than one chart plotted in one display, every chart has a position and a size
        chart "Lotka Volterra Time Series - Integration time step = 1 " type: series background: #white position: {0,0} size:{1,0.33} x_range: 1000 {       // one chart, of type 'serie', is named Lotka Volterra Time Series - Integration time step = 1, it shows quantities according to time, and the background is white
            data 'Number of preys' value: first(LotkaVolterra_agent where (each.integration_time_step = 1.0)).nb_prey color: #green ;           // number of preys in the case where the integration time step is 1 is plotted in green     
            data 'Number of predators' value: first(LotkaVolterra_agent where (each.integration_time_step = 1.0)).nb_predator color: #red ;     // number of predators in the case where the integration time step is 1 is plotted in red   
        }
        chart "Lotka Volterra Time Series - Integration time step = 0.1 " type: series background: #white position: {0,0.33} size:{1,0.33} x_range: 1000{
            data 'Number of preys' value: first(LotkaVolterra_agent where (each.integration_time_step = 0.1)).nb_prey color: #green ;               
            data 'Number of predators' value: first(LotkaVolterra_agent where (each.integration_time_step = 0.1)).nb_predator color: #red ;
        }
        chart "Lotka Volterra Time Series - Integration time step = 0.01 " type: series background: #white position: {0,0.66} size:{1,0.33}x_range: 1000{
            data 'Number of preys' value: first(LotkaVolterra_agent where (each.integration_time_step = 0.01)).nb_prey color: #green ;              
            data 'Number of predators' value: first(LotkaVolterra_agent where (each.integration_time_step = 0.01)).nb_predator color: #red ;
        }
    }
    display PhasePortrait {         
        chart "Lotka Volterra Phase Portrait - Integration time step = 1" type: xy background: #white position: {0,0} size:{1,0.33} {       // creation of a display to show phase portrait of the model, values are plotted at every step. Since there is more than one chart plotted in one display, every chart has a position and a size
        data 'Number of preys according to number of predators' value:{LV_agents[0].nb_prey, LV_agents[0].nb_predator} color: #black ;  // number of predators are plotted in black according to the number of preys in the case where the integration time step is 1       
        }
        chart "Lotka Volterra Phase Portrait - Integration time step = 0.1" type: xy background: #white position: {0,0.33} size:{1,0.33}{
        data 'Number of preys according to number of predators' value:{LV_agents[1].nb_prey, LV_agents[1].nb_predator} color: #black ;              
        }
        chart "Lotka Volterra Phase Portrait - Integration time step = 0.01" type: xy background: #white position: {0,0.66} size:{1,0.33} {
        data 'Number of preys according to number of predators' value:{LV_agents[2].nb_prey, LV_agents[2].nb_predator} color: #black ;              
        }
    }
}

}

AlexisDrogoul commented 6 years ago

Are these models committed ? Can we close this issue ?