#include <stdio.h>
#include <math.h>
#define x0 2.0
#define v0 1.0
#define t0 0.0
#define tf 100.0
#define w02 3.0
#define w0 (sqrt(w02))
#define xi 1.1
#define F0m 0
#define omega w0*sqrt(1-0.2)
#define dt 1e-2
#define N ((tf-t0)/dt)
double phi(double x, double t, double v)
{
return -w02*x + F0m*sin(omega*t) - 2*xi*v*w0;
}
double E(double x, double v)
{
return v*v/w02 + x*x;
}
int main(void)
{
double x = x0, xn,
t = t0, tn,
v = v0, vn;
FILE *h = fopen("./diferencial.dat","w+");
double E0 = E(x,v);
for (int i = 0; i < N; ++i)
{
tn = t + dt;
xn = x + vn*dt + 0.5*phi(x,t,v)*dt*dt;
vn = v + 0.5*(phi(xn,tn,v)+phi(x,t,v))*dt;
fprintf(h, "%d %f %f %f\n", i, t, x, E(x,v)/E0-1);
t = tn;
x = xn;
v = vn;
}
fclose(h);
return 0;
}
Comments