Ir al contenido principal

Nicolás' blog

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

Nicolás Guarín-Zapata

Comentarios

También disponible en:

English

Durante octubre (2017) estaré escribiendo un programa por día para algunos métodos numéricos famosos en Python y Julia. Esto está pensado como un ejercicio, no esperen que el código sea lo suficientemente bueno para usarse en la "vida real". Además, también debo mencionar que casi que no tengo experiencia con Julia, así que probablemente no escriba un Julia idiomático y se parezca más a Python.

Regla de Simpson

Hoy tenemos la regla de Simpson. Esta es otra de las fórmulas de Newton-Cotes, que se usa para la integración numérica. Newton-Cotes está relacionada con la interpolación de Lagrange con puntos equidistantes.

A continuación se presentan los códigos.

Python

from __future__ import division, print_function
import numpy as np


def simps(y, x=None, dx=1.0):
    n = len(y)
    if x is None:
        inte = np.sum(y[0:n-2:2] + 4*y[1:n-1:2] + y[2:n:2])*dx
    else:
        dx = x[1::2] - x[:-1:2]
        inte = np.dot(y[0:n-2:2] + 4*y[1:n-1:2] + y[2:n:2], dx)
    return inte/3


n = 21
x = np.linspace(0, 10, n)
y = np.sin(x)
print(simps(y, x=x))
print(1 - np.cos(10))

con resultados

1.8397296125
1.83907152908

Julia

function simps(y; x=nothing, dx=1.0)
    n = length(y)
    if x == nothing
        inte = sum(y[1:2:n-2] + 4*y[2:2:n-1] + y[3:2:n])*dx
    else
        dx = x[2:2:end] - x[1:2:end-1]
        inte = (y[1:2:n-2] + 4*y[2:2:n-1] + y[3:2:n])' * dx
    end
    return inte/3
end


n = 21
x = linspace(0, 10, n)
y = sin.(x)
println(simps(y, x=x))
println(1 - cos(10))

con resultados

1.8397296124951257
1.8390715290764525

Comparación Python/Julia

Respecto al número de líneas tenemos: 19 en Python y 17 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 simps(y, x=x)

con resultado

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

Para Julia:

@benchmark simps(y, x=x)

con resultado

BenchmarkTools.Trial:
  memory estimate:  1.23 KiB
  allocs estimate:  14
  --------------
  minimum time:     1.117 μs (0.00% GC)
  median time:      1.200 μs (0.00% GC)
  mean time:        1.404 μs (7.04% GC)
  maximum time:     222.286 μs (96.45% GC)
  --------------
  samples:          10000
  evals/sample:     10

En este caso, podemos decir que el código de Python es alrededor de 10 veces más lento que el de Julia.

Comentarios

Comments powered by Disqus