Skip to main content

Nicolás' blog

Numerical methods challenge: Day 13

Nicolás Guarín-Zapata

Comments

Also available in:

Español

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.

Cubic splines

Today we have cubic spline interpolation. Cubic splines are commonly used because it provides an approximation of a function with continuity up to second derivatives. For more details you can check this algorithm. Our main difference is that we will form the matrix and then solve the system. It is more common to directly solve the system of equation since it is a tridiagonal one.

Following are the codes.

Python

from __future__ import division, print_function
import numpy as np
import matplotlib.pyplot as plt


def spline_coeff(x, y):
    h = x[1:] - x[:-1]
    d_up = h.copy()
    d_up[0] = 0
    d_down = h.copy()
    d_down[-1] = 0
    d_cent = np.ones_like(x)
    d_cent[1:-1] = 2*(h[:-1] + h[1:])
    mat = np.diag(d_cent) + np.diag(d_up, 1) + np.diag(d_down, -1)
    alpha = np.zeros_like(x)
    alpha[1:-1] = 3/h[1:]*(y[2:] - y[1:-1]) - 3/h[:-1]*(y[1:-1] - y[:-2])
    c = np.linalg.solve(mat, alpha)
    b = np.zeros_like(x)
    d = np.zeros_like(x)
    b[:-1] = (y[1:] - y[:-1])/h - h/3*(c[1:] + 2*c[:-1])
    d[:-1] = (c[1:] - c[:-1])/(3*h)
    return y, b, c, d


n = 20
x = np.linspace(0, 2*np.pi, n)
y = np.sin(x)
a, b, c, d = spline_coeff(x, y)
for cont in range(n - 1):
    x_loc = np.linspace(x[cont], x[cont + 1], 100)
    y_loc = a[cont] + b[cont]*(x_loc - x[cont]) +\
            c[cont]*(x_loc - x[cont])**2 +\
            d[cont]*(x_loc - x[cont])**3
    plt.plot(x_loc, y_loc, "red")
plt.plot(x, y, "o")
plt.show()

Julia

using PyPlot


function spline_coeff(x, y)
    h = x[2:end] - x[1:end-1]
    d_up = copy(h)
    d_up[1] = 0.0
    d_down = copy(h)
    d_down[end-1] = 0.0
    d_cent = ones(x)
    d_cent[2:end-1] = 2*(h[1:end-1] + h[2:end])
    mat = diagm(d_cent) + diagm(d_up, 1) + diagm(d_down, -1)
    alpha = zeros(x)
    alpha[2:end-1] = 3./h[2:end].*(y[3:end] - y[2:end-1]) - 3./h[1:end-1].*(y[2:end-1] - y[1:end-2])
    c = mat \ alpha
    b = zeros(x)
    d = zeros(x)
    b[1:end-1] = (y[2:end] - y[1:end-1])./h - h./3.*(c[2:end] + 2*c[1:end-1])
    d[1:end-1] = (c[2:end] - c[1:end-1])./(3*h)
    return y, b, c, d
end


n = 21
x = collect(linspace(0, 2*pi, n))
y = sin.(x)
a, b, c, d = spline_coeff(x, y)
for cont = 1:n - 1
    x_loc = linspace(x[cont], x[cont + 1], 100)
    y_loc = a[cont] + b[cont]*(x_loc - x[cont]) +
            c[cont]*(x_loc - x[cont]).^2 +
            d[cont]*(x_loc - x[cont]).^3
    plot(x_loc, y_loc, "red")
end
plot(x, y, "o")
show()

In both cases the result is the plot below.

Spline interpolation.

Comparison Python/Julia

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

For Python:

%timeit a, b, c, d = spline_coeff(x, y)

with result

1000 loops, best of 3: 216 µs per loop

For Julia:

@benchmark a, b, c, d = spline_coeff(x, y)

with result

BenchmarkTools.Trial:
  memory estimate:  31.59 KiB
  allocs estimate:  52
  --------------
  minimum time:     18.024 μs (0.00% GC)
  median time:      26.401 μs (0.00% GC)
  mean time:        44.035 μs (3.94% GC)
  maximum time:     9.833 ms (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1

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

Comments

Comments powered by Disqus