I wrote a simple code implementing a midpoint 2nd order Runge-Kutta method for a damped driven pendulum in different languages to compare run times. While the results were exactly the same for C, Fortran and Python (Fortran and C tested up to t = 10⁷) the results for Julia drifted away from these after a rather short time.
Since the problem is chaotic with the given parameters, any slight difference the initial conditions will grow with time; however, all languages tested except for Julia apparently were able to perform exactly the same calculations.
Can someone explain what is happening? (Equations, Julia and C codes below image)
\begin{cases}\frac{d\theta}{dt} = v \\
\frac{dv}{dt} = -\gamma\, v - \sin(\theta) + F \cos(\omega\, t)
\end{cases}
x_0 = 1.0,\; v_0 = 0.0
\omega = 0.7, \; \gamma = 0.05, \; F = 0.6
Julia code:
using Printf
X0 = 1.0
VX0 = 0.0
OMEGA = 0.70
GAMMA = 0.05
FF = 0.6
TMAX = 500
DT = 0.01
DT_DIV_2 = 0.5*DT
function f(t, x, v)
return -GAMMA*v - sin(x) + FF*cos(OMEGA*t)
end
out1 = @sprintf("dados_rk2_F%.2f_julia.txt", FF)
fout1 = open(out1, "w")
i = 0
x = X0
v = VX0
t = 0.0
while (t<=TMAX)
k1x = v
k1v = f(t,x,v)
tt = t + DT_DIV_2
xx = x + k1x * DT_DIV_2
vv = v + k1v * DT_DIV_2
k2x = vv
k2v = f(tt, xx, vv)
global x += k2x * DT
global v += k2v * DT
global i += 1
global t = i*DT
@printf(fout1,"%12.6f %26.20f %26.20f\n", t, x, v)
end
close(fout1)
C language code:
#include <stdio.h>
#include <math.h>
const double X0 = 1.0;
const double VX0 = 0.0;
const double OMEGA = 0.70;
const double GAMMA = 0.05;
const double FF = 0.60;
const double TMAX = 500;
const double DT = 0.01;
double f(double t, double x, double v)
{
return -GAMMA*v - sin(x) + FF*cos(OMEGA*t);
}
int main(int argc, char *argv[])
{
double x, v, t;
double xx, vv, tt, k1x, k2x, k1v, k2v;
double DT_DIV_2 = 0.5*DT;
char out1[200];
FILE *fout1;
int i = 0;
sprintf(out1,"dados_rk2_F%.2f.txt",FF);
fout1=fopen(out1,"w");
x = X0; v = VX0; t = 0.0;
while (t <= TMAX)
{
k1x = v;
k1v = f(t,x,v);
tt = t + DT_DIV_2;
xx = x + k1x * DT_DIV_2;
vv = v + k1v * DT_DIV_2;
k2x = vv;
k2v = f(tt, xx, vv);
x += k2x * DT;
v += k2v * DT;
i++;
t = i*DT;
fprintf(fout1,"%12.6f %26.20f %26.20f\n", t, x, v);
}
fclose(fout1);
return 0;
}
