Reto de métodos numéricos: Día 14

Python

from __future__ import division, print_function
import numpy as np


def trapz(y, x=None, dx=1.0):
    if x is None:
        inte = 0.5*dx*(2*np.sum(y[1:-1]) + y[0] + y[-1])
    else:
        dx = x[1:] - x[:-1]
        inte = 0.5*np.dot(y[:-1] + y[1:], dx)
    return inte


dx = 0.001
x = np.arange(0, 10, dx)
y = np.sin(x)
print(trapz(y, dx=dx))
print(1 - np.cos(10))

con resultado

1.83961497726
1.83907152908

Julia

function trapz(y; x=nothing, dx=1.0)
    if x == nothing
        inte = 0.5*dx*(2*sum(y[2:end-1]) + y[1] + y[end])
    else
        dx = x[2:end] - x[1:end-1]
        inte = 0.5* (y[1:end-1] + y[2:end])' * dx
    end
    return inte
end


dx = 0.001
x = 0:dx:10
y = sin.(x)
println(trapz(y, dx=dx))
println(1 - cos(10))

con resultado

1.8390713758204895
1.8390715290764525

Comparación Python/Julia

Respecto al número de líneas tenemos: 18 en Python y 16 en Julia. La comparación en tiempo de ejecución se realizó con el comando mágico de IPython %timeit y con @benchmark en Julia.

Para Python:

%timeit trapz(y, dx=dx)

con resultado

100000 loops, best of 3: 16.9 µs per loop

Para Julia:

@benchmark trapz(y, dx=dx)

con resultado

BenchmarkTools.Trial:
  memory estimate:  78.31 KiB
  allocs estimate:  4
  --------------
  minimum time:     13.080 μs (0.00% GC)
  median time:      16.333 μs (0.00% GC)
  mean time:        20.099 μs (12.66% GC)
  maximum time:     963.732 μs (90.60% GC)
  --------------
  samples:          10000
  evals/sample:     1

En este caso, podemos decir que el código de Python es tan rápido como el de Julia.