# Numerical methods challenge: Day 15

During October (2017) I will write a program per day for some well-known numerical methods in both Python and Julia. It is intended to be an exercise then don't expect the code to be good enough for real use. Also, I should mention that I have almost no experience with Julia, so it probably won't be idiomatic Julia but more Python-like Julia.

## Simpson's rule

Today we have the Simpson's rule. This is another Newton-Cotes formulas, used for numerical integration. Newton-Cotes is related to Lagrange interpolation with equidistant points.

Following are the codes.

### 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))

with result

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))

1.8397296124951257 1.8390715290764525

### Comparison Python/Julia

Regarding number of lines we have: 19 in Python and 17 in Julia. The comparison
in execution time is done with `%timeit` magic command in IPython and
`@benchmark` in Julia.

For Python:

%timeit simps(y, x=x)

with result

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

For Julia:

@benchmark simps(y, x=x)

with result

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

In this case, we can say that the Python code is 10 times slower than Julia.

## Comments

Comments powered by Disqus