# Numerical difference between Julia and (C, Fortran, Python) for a chaotic damped driven pendulum

**URL:** <https://discourse.julialang.org/t/numerical-difference-between-julia-and-c-fortran-python-for-a-chaotic-damped-driven-pendulum/62195>\
**Category:** Modelling & Simulations\
**Tags:** numerics, precision\
**Created:** [June 1, 2021, 2:42pm UTC](https://discourse.julialang.org/t/numerical-difference-between-julia-and-c-fortran-python-for-a-chaotic-damped-driven-pendulum/62195 "2021-06-01T14:42:33Z")\
**Posts on this page:** 18\
**Page:** 1

<div class="post-metadata">

**Author:** ![mendeli](https://avatars.discourse-cdn.com/v4/letter/m/6de8d8/32.png) [@mendeli](https://discourse.julialang.org/u/mendeli)\
**Post date:** [June 1, 2021, 2:42pm UTC](https://discourse.julialang.org/t/numerical-difference-between-julia-and-c-fortran-python-for-a-chaotic-damped-driven-pendulum/62195/1 "2021-06-01T14:42:33Z")

</div>

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)

 ![Screenshot from 2021-06-01 09-27-15](https://global.discourse-cdn.com/julialang/original/3X/6/6/6657a3d1769045dd6db9a3985bbaaa6d5cc0cb12.png)

\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:

```julia
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:

```julia
#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;
}

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 1, 2021, 2:43pm UTC](https://discourse.julialang.org/t/numerical-difference-between-julia-and-c-fortran-python-for-a-chaotic-damped-driven-pendulum/62195/2 "2021-06-01T14:43:52Z")

</div>

> [@mendeli](#):
>
> I wrote a simple code implementing a midpoint 2nd order Runge-Kutta method for a damped driven pendulum in different languages to compared 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.

It’s chaotic, so the simplest things can cause this. My guess is that the first 3 use the same standard math library definition of `sin` and `cos`, while Julia’s is different in 1ulp or so.

---

<div class="post-metadata">

**Author:** ![mendeli](https://avatars.discourse-cdn.com/v4/letter/m/6de8d8/32.png) [@mendeli](https://discourse.julialang.org/u/mendeli)\
**Post date:** [June 1, 2021, 2:56pm UTC](https://discourse.julialang.org/t/numerical-difference-between-julia-and-c-fortran-python-for-a-chaotic-damped-driven-pendulum/62195/3 "2021-06-01T14:56:55Z")

</div>

Thanks for the heads up. I naively expected that all implementations for the transcendental functions should work exactly the same with double precision. If anyone is curious:

> **[Unit in the last place](https://en.wikipedia.org/wiki/Unit_in_the_last_place)**
>
> In computer science and numerical analysis, unit in the last place or unit of least precision (ulp) is the spacing between two consecutive floating-point numbers, i.e., the value the least significant digit (rightmost digit) represents if it is 1. It is used as a measure of accuracy in numeric calculations.
> One definition is: In radix 
>   
>     
>       
> b
>       
>     
> {\\displaystyle b}
>   
> with precision 
>   
>     
>       
> p
>       
>     
> {\\displaystyle p}
>   
> , if 
>   
>     
> ...

“Reputable numeric libraries compute the basic transcendental functions to between 0.5 and about 1 ULP. Only a few libraries compute them within 0.5 ULP, this problem being complex due to the Table-maker’s dilemma.”

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [June 1, 2021, 2:57pm UTC](https://discourse.julialang.org/t/numerical-difference-between-julia-and-c-fortran-python-for-a-chaotic-damped-driven-pendulum/62195/4 "2021-06-01T14:57:42Z")

</div>

You can make some rough estimate of the divergence due to the chaotic dynamics based on the Lyapunov exponent. Here is a simple script that calculates the maximum Lyapunov exponent:

```julia
using DynamicalSystems

function pendulum_rule(u, p, t)
    ω, γ, F = p
    θ, v = u
    dθ = v
    dv = -γ*v - sin(θ) + F*cos(ω*t)
    return SVector(dθ, dv)
end

p0 = [0.7, 0.05, 0.6]
u0 = [1.0 + 1e-3rand(), 0.0]

ds = ContinuousDynamicalSystem(pendulum_rule, u0, p0)

λ = lyapunov(ds, 100000.0)

e = 1e-12 # error expected 

Δt = log(1/e)/λ

```

```julia
195.5...

```

Using this I get a prediction horizon of about 200 time units, which seems to go well with your plot. Of course, this is not precisely accurate, as the roundoff error due to cosine differences is more similar to dynamic noise constantly added to the system, but can give you a good estimate of how fast you would expect divergence.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [June 1, 2021, 3:14pm UTC](https://discourse.julialang.org/t/numerical-difference-between-julia-and-c-fortran-python-for-a-chaotic-damped-driven-pendulum/62195/5 "2021-06-01T15:14:05Z")

</div>

You can try the system `sin` and `cos` functions instead of Julia’s:

```julia
csin(x::Float64) = ccall(:sin, Float64, (Float64,), x)
csin(x::Float32) = ccall(:sinf, Float32, (Float32,), x)
ccos(x::Float64) = ccall(:cos, Float64, (Float64,), x)
ccos(x::Float32) = ccall(:cosf, Float32, (Float32,), x)

```

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [June 1, 2021, 3:20pm UTC](https://discourse.julialang.org/t/numerical-difference-between-julia-and-c-fortran-python-for-a-chaotic-damped-driven-pendulum/62195/6 "2021-06-01T15:20:03Z")

</div>

> [@mendeli](#):
>
> to compare run times.

Before you get disapointed, your code as is will be very slow in Julia. Put all inside a function and benchmark it with BenchmarkTools to actually get its peformance. Something like this (be aware that in all these benchmarks the time will be dominated by printing as they are):

```julia
using Printf, BenchmarkTools

function pendulum(;print=false)

  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

  f(t, x, v) = -GAMMA*v - sin(x) + FF*cos(OMEGA*t)

  if print
    out1 = @sprintf("dados_rk2_F%.2f_julia.txt", FF)
    fout1 = open(out1, "w")
  end

  i = 0
  x = X0 
  t = 0.0
  v = VX0
  
  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 += 1
  	t = i*DT
  		
     if print 
       @printf(fout1,"%12.6f %26.20f %26.20f\n", t, x, v)
     end
  end
  if print
    close(fout1)
  end
	
end

pendulum(print=true)

@benchmark pendulum(print=false)

```

You will get:

```julia
julia> include("./pend.jl")
BenchmarkTools.Trial:
  memory estimate: 0 bytes
  allocs estimate: 0
  --------------
  minimum time: 1.795 ms (0.00% GC)
  median time: 1.868 ms (0.00% GC)
  mean time: 1.875 ms (0.00% GC)
  maximum time: 2.643 ms (0.00% GC)
  --------------
  samples: 2662
  evals/sample: 1

```

---

<div class="post-metadata">

**Author:** ![mendeli](https://avatars.discourse-cdn.com/v4/letter/m/6de8d8/32.png) [@mendeli](https://discourse.julialang.org/u/mendeli)\
**Post date:** [June 1, 2021, 5:06pm UTC](https://discourse.julialang.org/t/numerical-difference-between-julia-and-c-fortran-python-for-a-chaotic-damped-driven-pendulum/62195/7 "2021-06-01T17:06:04Z")

</div>

I tried using the Float64 versions in my program but I still got the same result as the original Julia program.

---

<div class="post-metadata">

**Author:** ![mendeli](https://avatars.discourse-cdn.com/v4/letter/m/6de8d8/32.png) [@mendeli](https://discourse.julialang.org/u/mendeli)\
**Post date:** [June 1, 2021, 5:24pm UTC](https://discourse.julialang.org/t/numerical-difference-between-julia-and-c-fortran-python-for-a-chaotic-damped-driven-pendulum/62195/8 "2021-06-01T17:24:14Z")

</div>

Thanks for the reply. Even if I switch off printing, the Julia code is still quite slower than C. Is there anything that can be done to improve performance, like Numba for Python?

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [June 1, 2021, 5:26pm UTC](https://discourse.julialang.org/t/numerical-difference-between-julia-and-c-fortran-python-for-a-chaotic-damped-driven-pendulum/62195/9 "2021-06-01T17:26:36Z")

</div>

There’s no need for anything like Numba–most likely you just need to follow @lmiq’s suggestion to put your code into a function (see [Performance Tips · The Julia Language](https://docs.julialang.org/en/v1/manual/performance-tips/#Avoid-global-variables) ). If that’s still not enough, definitely post your updated code and its performance here, and I’m sure we can find some more concrete suggestions.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [June 1, 2021, 5:30pm UTC](https://discourse.julialang.org/t/numerical-difference-between-julia-and-c-fortran-python-for-a-chaotic-damped-driven-pendulum/62195/10 "2021-06-01T17:30:25Z")

</div>

> [@mendeli](#):
>
> Even if I switch off printing, the Julia code is still quite slower than C.

As I posted above, you need to put your code inside a function (see the code I posted there). Also, Julia compiles the code the first time it is run, so for very quick tests like this you might be measuring the compilation time as an important part of the total time. That is why using the benchmark tools is important there.

With both codes printing (which is almost all the time), here the Julia code (I posted) takes 49ms and the C code (measured with system’s time, thus probably not a good measure either), takes 85ms. But with good performance measures I would expect both to be similar here.

Ps: If you want a simpler but reasonably fair comparison, take the code I posted there, remove the last line (`@benchmark...`) and change the number of steps to something much bigger (1 million). Probably removing printing is a good idea also:

```julia
%gcc -o pend pend.c -lm -O3
% time ./pend
real 0m6,144s
user 0m6,139s
sys 0m0,004s

% time julia ./pend.jl
real 0m4,854s
user 0m5,167s
sys 0m0,628s

```

With these codes:

> **Julia**
>
> ```julia
> using Printf, BenchmarkTools
> 
> function pendulum(;print=false,TMAX=500)
> 
> X0 = 1.0		
> VX0 = 0.0		
>   
> OMEGA = 0.70
> GAMMA = 0.05
> FF = 0.6 
>   
> DT = 0.01  
> DT_DIV_2 = 0.5*DT
> 
> f(t, x, v) = -GAMMA*v - sin(x) + FF*cos(OMEGA*t)
> 
> if print
> out1 = @sprintf("dados_rk2_F%.2f_julia.txt", FF)
> fout1 = open(out1, "w")
> end
> i = 0
> x = X0 
> t = 0.0
> v = VX0
>   
> 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 += 1
> t = i*DT
>      	
> if print 
> @printf(fout1,"%12.6f %26.20f %26.20f\n", t, x, v)
> end
> end
> if print
> close(fout1)
> end 
> return x
> 
> end
> 
> pendulum(print=false,TMAX=1_000_000)
> 
> ```

> **C**
>
> ```julia
> #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 = 1000000;
> 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 x;
> }
> 
> ```

(note that I had to return `x` from both functions, otherwise the compiler with `-O3` just figures out that it does not need to do anything if not printing the trajectory)

---

<div class="post-metadata">

**Author:** ![mendeli](https://avatars.discourse-cdn.com/v4/letter/m/6de8d8/32.png) [@mendeli](https://discourse.julialang.org/u/mendeli)\
**Post date:** [June 1, 2021, 5:36pm UTC](https://discourse.julialang.org/t/numerical-difference-between-julia-and-c-fortran-python-for-a-chaotic-damped-driven-pendulum/62195/11 "2021-06-01T17:36:41Z")

</div>

Thank you @rdeits and @lmiq! This was my first attempt at using Julia.

---

<div class="post-metadata">

**Author:** ![mendeli](https://avatars.discourse-cdn.com/v4/letter/m/6de8d8/32.png) [@mendeli](https://discourse.julialang.org/u/mendeli)\
**Post date:** [June 1, 2021, 5:43pm UTC](https://discourse.julialang.org/t/numerical-difference-between-julia-and-c-fortran-python-for-a-chaotic-damped-driven-pendulum/62195/12 "2021-06-01T17:43:36Z")

</div>

That was very helpful! Thank you very much.

---

<div class="post-metadata">

**Author:** ![mendeli](https://avatars.discourse-cdn.com/v4/letter/m/6de8d8/32.png) [@mendeli](https://discourse.julialang.org/u/mendeli)\
**Post date:** [June 1, 2021, 5:58pm UTC](https://discourse.julialang.org/t/numerical-difference-between-julia-and-c-fortran-python-for-a-chaotic-damped-driven-pendulum/62195/13 "2021-06-01T17:58:47Z")

</div>

I’m impressed! Even printing every 100 steps it is running faster than C.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [June 1, 2021, 6:05pm UTC](https://discourse.julialang.org/t/numerical-difference-between-julia-and-c-fortran-python-for-a-chaotic-damped-driven-pendulum/62195/14 "2021-06-01T18:05:16Z")

</div>

On neither side you should expect one code to be much faster than the other. If one or other thing is happening, probably there is something to be improved on the other side. I am saying this so that expectations are reasonable: to get performance in Julia you need to get used to some few coding strategies that are inherent to it. And with that you will get code that is of the order of C-fast.

From that point on one or other code can get faster with more specialized optimizations, tuning compiler flags, etc.

Juila is nice because the experience of coding that is much smoother in general (once one gets used to it, as in all cases).

---

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [June 1, 2021, 6:40pm UTC](https://discourse.julialang.org/t/numerical-difference-between-julia-and-c-fortran-python-for-a-chaotic-damped-driven-pendulum/62195/15 "2021-06-01T18:40:49Z")

</div>

At least on my machine, gcc doesn’t want to use AVX, whereas Julia seemingly does.

I suspect this vectorization is what makes the difference here.

> **C asm**
>
> ```nohighlight
> 0000000000001420 <f>:                                                                   
> 1420: f3 0f 1e fa endbr64                                         
> 1424: 48 83 ec 28 sub $0x28,%rsp                               
> 1428: f2 0f 11 44 24 18 movsd %xmm0,0x18(%rsp)                         
> 142e: 66 0f 28 c1 movapd %xmm1,%xmm0                              
> 1432: f2 0f 11 54 24 10 movsd %xmm2,0x10(%rsp)                         
> 1438: e8 b3 fc ff ff callq 10f0 <sin@plt>                           
> 143d: f2 0f 10 5c 24 18 movsd 0x18(%rsp),%xmm3                         
> 1443: f2 0f 59 1d 25 0c 00 mulsd 0xc25(%rip),%xmm3 # 2070 <X0+0x8> 
> 144a: 00                                                                      
> 144b: f2 0f 11 44 24 08 movsd %xmm0,0x8(%rsp)                          
> 1451: 66 0f 28 c3 movapd %xmm3,%xmm0                              
> 1455: e8 76 fc ff ff callq 10d0 <cos@plt>                           
> 145a: f2 0f 59 05 1e 0c 00 mulsd 0xc1e(%rip),%xmm0 # 2080 <X0+0x18>
> 1461: 00                                                                      
> 1462: f2 0f 10 54 24 10 movsd 0x10(%rsp),%xmm2                         
> 1468: f2 0f 59 15 08 0c 00 mulsd 0xc08(%rip),%xmm2 # 2078 <X0+0x10>
> 146f: 00                                                                      
> 1470: f2 0f 5c 54 24 08 subsd 0x8(%rsp),%xmm2                          
> 1476: 48 83 c4 28 add $0x28,%rsp                               
> 147a: f2 0f 58 c2 addsd %xmm2,%xmm0                              
> 147e: c3 retq                                            
> 147f: 90 nop                                             
> 
> ```

vs.

> **Julia asm**
>
> ```julia
> julia> @code_native acceleration(0.0, 1.0, 0.0, 0.05, 0.6, 0.7)     
> .text                                                       
> ; ┌ @ REPL[42]:1 within `acceleration'                              
> pushq %rbx                                                
> subq $32, %rsp                                           
> vmovsd %xmm5, 16(%rsp)                                     
> vmovsd %xmm4, 24(%rsp)                                     
> vmovsd %xmm0, 8(%rsp)                                      
> movabsq $.rodata.cst16, %rax                                
> ; │┌ @ float.jl:383 within `-'                                      
> vxorpd (%rax), %xmm3, %xmm0                                
> ; │└                                                                
> ; │┌ @ float.jl:395 within `*'                                      
> vmulsd %xmm2, %xmm0, %xmm0                                 
> vmovsd %xmm0, (%rsp)                                       
> movabsq $cos, %rbx                                          
> ; │└                                                                
> ; │┌ @ REPL[38]:1 within `csin'                                     
> leaq 17120(%rbx), %rax                                   
> vmovapd %xmm1, %xmm0                                        
> callq *%rax                                               
> vmovsd (%rsp), %xmm1 # xmm1 = mem[0],zero
> ; │└                                                                
> ; │┌ @ float.jl:392 within `-'                                      
> vsubsd %xmm0, %xmm1, %xmm0                                 
> vmovsd %xmm0, (%rsp)                                       
> vmovsd 8(%rsp), %xmm0 # xmm0 = mem[0],zero
> ; │└                                                                
> ; │┌ @ float.jl:395 within `*'                                      
> vmulsd 16(%rsp), %xmm0, %xmm0                              
> ; │└                                                                
> ; │┌ @ REPL[40]:1 within `ccos'                                     
> callq *%rbx                                               
> ; │└                                                                
> ; │┌ @ float.jl:395 within `*'                                      
> vmulsd 24(%rsp), %xmm0, %xmm0                              
> ; │└                                                                
> ; │┌ @ float.jl:389 within `+'                                      
> vaddsd (%rsp), %xmm0, %xmm0                                
> ; │└                                                                
> addq $32, %rsp                                           
> popq %rbx                                                
> retq                                                        
> nopw %cs:(%rax,%rax)                                     
> ; └                                                                 
> 
> ```

no, the calls are not inlined in julia (not with `ccos` and not with native `sin`), though I did rename them for my understanding - `acceleration` is your `f`:

> **callq are in here**
>
> ````julia
> L1340:                                                              
> vmovsd %xmm0, 32(%rsp)                                     
> movabsq $sin, %rbx                                          
> ; │ @ chaotic.jl:18 within `pendulum'                               
> ; │┌ @ chaotic.jl:13 within `acceleration'                          
> callq *%rbx                                               
> vmovsd %xmm0, 120(%rsp)                                    
> vmovsd 24(%rsp), %xmm0 # xmm0 = mem[0],zero
> ; ││┌ @ float.jl:395 within `*'                                     
> vmulsd 80(%rsp), %xmm0, %xmm0                              
> movabsq $cos, %r15                                          
> ; ││└                                                               
> callq *%r15                                               
> vmovsd %xmm0, 112(%rsp)                                    
> vmovsd 72(%rsp), %xmm0 # xmm0 = mem[0],zero
> ; │└ ```
> 
> ````

> **Julia Code**
>
> ```julia
> function pendulum(print=false, x=1.0, v=0.0, dt=0.01, tmax=500.0, Ω=0.70, Γ=0.05, FF=0.6)
> if print
> out1 = @sprintf("2_dados_rk2_F%.2f_julia.txt", FF)
> fout1 = open(out1, "w")
> end
> acceleration(t, x, v) = -Γ * v - sin(x) + FF * cos(Ω * t)
> 
> dt_half = dt / 2.0
> 
> for t in range(0.0, tmax; step=dt)
> acc_start = acceleration(t, x, v)
>         
> t_middle = t + dt_half
> x_middle = x + v * dt_half
> v_middle = v + acc_start * dt_half	          
> acc_middle = acceleration(t_middle, x_middle, v_middle)                  
>         
> x += v_middle * dt	
> v += acc_middle * dt
>   
> if print
> @printf(fout1,"%12.6f %26.20f %26.20f\n", t+dt, x, v)
> end
> end
> if print
> close(fout1)
> end
> end
> 
> ```

C code was the one by @lmiq, compiled with `gcc -o pendulum pendulum.c -lm -O3`.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [June 1, 2021, 8:45pm UTC](https://discourse.julialang.org/t/numerical-difference-between-julia-and-c-fortran-python-for-a-chaotic-damped-driven-pendulum/62195/16 "2021-06-01T20:45:14Z")

</div>

Now the really cool thing is this. You can with very minor modifications to your function make simulations of the pendulums of any dimension, defined by the input type of the initial point and velocities. Now do that with C 🙂

 ![pend](https://global.discourse-cdn.com/julialang/original/3X/4/3/43e9a9a1851cb9ae713b9360b95c0f502fd89c27.png)

There are a few annotations on where I had to change the code (besides the plotting, of course):

```julia
using Printf, StaticArrays

function pendulum(X0,VX0;printout=false,TMAX=500) # X0 and VX0 are input parameters

  OMEGA = 0.70
  GAMMA = 0.05
  FF = 0.6 
  
  DT = 0.01  
  DT_DIV_2 = 0.5*DT

  f(t, x, v) = -GAMMA*v - sin(x) + FF*cos(OMEGA*t)

  if printout
    out1 = "dados_julia.txt"
    fout1 = open(out1, "w")
  end
  i = 0
  x = X0 
  t = 0.0
  v = VX0
  
  while t <= TMAX

    k1x = v
    k1v = f.(t,x,v) # Added a dot here
    	
    tt = t + DT_DIV_2
    xx = x + k1x * DT_DIV_2
    vv = v + k1v * DT_DIV_2	          
    	
    k2x = vv
    k2v = f.(tt, xx, vv) # added a dot here                  
    	  
    x += k2x * DT	
    v += k2v * DT
    	
    i += 1
    t = i*DT
     	
    if printout && i%5 == 0 # print generic sizes 
      println(fout1,"$t ",["$xi " for xi in x]...)
    end
  end
  if printout
    close(fout1)
  end
  return x

end

using DelimitedFiles, Plots, Plots.Measures

plot(layout=(1,3))
default(label="")

# 1D pendulum
x0 = 1.0
vx0 = rand()
pendulum(x0,vx0,printout=true,TMAX=1_000)
x1 = readdlm("dados_julia.txt")
plot!(x1[:,1],x1[:,2],subplot=1,title="1D pendulum")
plot!(xlabel="time",ylabel="x",subplot=1)

# 2D pendulum
n = 2
x0 = ones(SVector{n,Float64})
vx0 = rand(SVector{n,Float64})
pendulum(x0,vx0,printout=true,TMAX=1_000)
x2 = readdlm("dados_julia.txt")
plot!(x2[:,1],x2[:,2],x2[:,3],subplot=2,title="2D pendulum")
plot!(xlabel="time",ylabel="x1",zlabel="x2",subplot=2)

# 3D pendulum
n = 3
x0 = ones(SVector{n,Float64})
vx0 = rand(SVector{n,Float64})
pendulum(x0,vx0,printout=true,TMAX=1_000)
x3 = readdlm("dados_julia.txt")
plot!(x3[:,2],x3[:,3],x3[:,4],subplot=3,title="3D pendulum")
plot!(xlabel="x1",ylabel="x2",zlabel="x3",subplot=3)

plot!(size=(1200,500),margin=5mm)
savefig("./pend.png")

```

---

<div class="post-metadata">

**Author:** ![mendeli](https://avatars.discourse-cdn.com/v4/letter/m/6de8d8/32.png) [@mendeli](https://discourse.julialang.org/u/mendeli)\
**Post date:** [June 2, 2021, 11:45am UTC](https://discourse.julialang.org/t/numerical-difference-between-julia-and-c-fortran-python-for-a-chaotic-damped-driven-pendulum/62195/18 "2021-06-02T11:45:40Z")

</div>

That’s really cool! Thank for the nice examples. I guess I’ll be learning Julia.

---

<div class="post-metadata">

**Author:** ![murrayE](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/murraye/32/24428_2.png) [@murrayE](https://discourse.julialang.org/u/murrayE)\
**Post date:** [June 8, 2021, 12:22am UTC](https://discourse.julialang.org/t/numerical-difference-between-julia-and-c-fortran-python-for-a-chaotic-damped-driven-pendulum/62195/19 "2021-06-08T00:22:17Z")

</div>

For what it’s worth, you may wish to compare the results also with those from Mathematica, which is presumably using state-of-the-art numerical methods.

Code:

 ![pendulum_code](https://global.discourse-cdn.com/julialang/original/3X/6/f/6fa7cac606d25f8be4074f4e16d3dd901c99e454.png)

Plot:

 ![pendulum_plot](https://global.discourse-cdn.com/julialang/original/3X/e/2/e273aacf211e22df987894adb643206d99905f49.png)