Hello guys, im a beginner programer and I need help with my simulation of the Lotka-Volterra model. I have to implement Euler and Verlet algorithm and I dont know if it is well done. In the model there is constant, V, my problem is that V is not conserved in time and I dont know why. I don't know what problems my code has, nor do I know if the results it gives me so far are correct. I need help, I leave the code here, thank you very much.
Code:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define dt 0.0625
#define MAX 400
#define EULER
void feuler(double xn,double yn,double *xnext,double *ynext);
void fverlet(double xn,double yn,double *xnext,double *ynext);
double Cxx=1.0;
double Cxy=0.2;
double Cyy=0.2;
double Cyx=0.01;
int main()
{
double x;
double y;
double t;
double xaux,yaux;
double V,Vaux;
FILE *f;
x=20;
y=6;
t=0;
f=fopen("PyD.txt","w");
if (f==NULL)
{
printf("ERROR");
}
else{
Vaux=Cyx*x-Cyy*log(x)+Cxy*y-Cxx*log(y);
fprintf(f," t\t\t x\t\t y\t\t V\t\t dx_dt\t dy_dt\n%lf\t%lf\t%lf\t%lf\t%lf\t %lf\n",t,x,y,Vaux,Cxx*x-Cxy*x*y,Cyx*x*y-Cyy*y);
#ifdef EULER
for(t=0;t<=MAX;t+=dt) // Este bucle representa el "tiempo" transcurrido.
{
feuler(x,y,&xaux,&yaux);
x=xaux;
y=yaux;
V=Cyx*x-Cyy*log(x)+Cxy*y-Cxx*log(y);
fprintf(f,"%lf\t%lf\t%lf\t%lf\n",t,x,y,V);
}
#else
for(t=0;t<MAX;t+=dt)
{
fverlet(x,y,&xaux,&yaux);
x=xaux;
y=yaux;
V=Cyx*x-Cyy*log(x)+Cxy*y-Cxx*log(y):
fprintf(f,"%lf\t%lf\t%lf\t%lf\n",t,x,y,V);
}
#endif
V=Cyx*x-Cyy*log(x)+Cxy*y-Cxx*log(y);
printf("V_i=%lf\nV_f=%lf\n",Vaux,V);
}
fclose(f);
return 0;
}
void feuler(double xn,double yn,double *xnext,double *ynext)
{
double dx_dt,dy_dt;
dx_dt=Cxx*xn-Cxy*xn*yn;
dy_dt=Cyx*xn*yn-Cyy*yn;
*xnext=xn+dx_dt*dt;
*ynext=yn+dy_dt*dt;
}
void fverlet(double xn,double yn,double *xnext,double *ynext)
{
double dx_dt,dy_dt;
double dx_dt2,dy_dt2;
double xmed,ymed;
dx_dt=Cxx*xn-Cxy*xn*yn;
dy_dt=Cyx*xn*yn-Cyy*yn;
dx_dt2=Cxx*dx_dt-Cxy+(dy_dt*xn+yn*dx_dt);
dy_dt2=Cyx*(dy_dt*xn+yn*dx_dt)-Cyy*dy_dt;
*xnext=xn+dx_dt*dt+dx_dt2*dt*dt*0.5;
*ynext=yn+dy_dt*dt+dy_dt2*dt*dt*0.5;
return;
}
Last edited: