Skip to main content

Numerical methods challenge: Day 2

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.

Regula falsi

The second method to be considered is the false position method, or regula falsi. This method is used to solve the equation \(f(x) = 0\) for \(x\) real, and \(f\) continuous. It starts with an interval \([a,b]\), where \(f(a)\) and \(f(b)\) should have opposite signs. The method is similar to bisection method but instead of halving the original interval it takes as a new point of the interval the intercept of the line that connect the function evaluated at the two extreme points. Then, the new point is computed from

\begin{equation*} c = \frac{a f(b) - b f(a)}{f(b) - f(a)} \end{equation*}

As in bisection method, we keep the interval where the solution appears, i.e., the sign of the function changes.

We will use the function \(f(x) = \cos(x) - x\) to test the codes, and the initial interval is \([0.5, \pi/4]\).

Following are the codes.

Python

from __future__ import division, print_function
from numpy import abs, cos, pi

def regula_falsi(fun, a, b, niter=50, ftol=1e-12, verbose=False):
    if fun(a) * fun(b) > 0:
        c = None
        msg = "The function should have a sign change in the interval."
    else:
        for cont in range(niter):
            qa = fun(a)
            qb = fun(b)
            c = (a*qb - b*qa)/(qb - qa)
            qc = fun(c)
            if verbose:
                print("n: {}, c: {}".format(cont, c))
            msg = "Maximum number of iterations reached."
            if abs(qc) < ftol:
                msg = "Root found with desired accuracy."
                break
            elif qa * qc < 0:
                b = c
            elif qb * qc < 0:
                a = c
    return c, msg

def fun(x):
    return cos(x) - x

print(regula_falsi(fun, 0.5, 0.25*pi))

Julia

function regula_falsi(fun, a, b, niter=50, ftol=1e-12, verbose=false)
    if fun(a) * fun(b) > 0
        c = nothing
        msg = "The function should have a sign change in the interval."
    else
        for cont = 1:niter
            qa = fun(a)
            qb = fun(b)
            c = (a*qb - b*qa)/(qb - qa)
            qc = fun(c)
            if verbose
                println("n: $(cont), c: $(c)")
            end
            if abs(fun(c)) < ftol
                msg = "Root found with desired accuracy."
                break
            elseif qa * qc < 0
                b = c
            elseif qb * qc < 0
                a = c
            end
            msg = "Maximum number of iterations reached."
        end
    end
    return c, msg
end

function fun(x)
    return cos(x) - x
end

println(regula_falsi(fun, 0.5, 0.25*pi))

Comparison

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

For Python:

%timeit regula_falsi(fun, 0.5, 0.25*pi)

with result

10000 loops, best of 3: 35.1 µs per loop

For Julia:

@benchmark regula_falsi(fun, 0.5, 0.25*pi)

with result

BenchmarkTools.Trial:
  memory estimate:  48 bytes
  allocs estimate:  2
  --------------
  minimum time:     449.495 ns (0.00% GC)
  median time:      464.371 ns (0.00% GC)
  mean time:        493.785 ns (0.52% GC)
  maximum time:     9.723 μs (92.54% GC)
  --------------
  samples:          10000
  evals/sample:     198

Numerical methods challenge

For the next month 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.

Day 1: Bisection method

The first method to be considered is the bisection method. This method is used to solve the equation \(f(x) = 0\) for \(x\) real, and \(f\) continuous. It starts with an interval \([a,b]\), where \(f(a)\) and \(f(b)\) should have opposite signs. The method then proceeds by halving the interval and selecting the one where the solution appears, i.e., the sign of the function changes.

We will use the function \(f(x) = \cos(x) - x^2\) to test the codes, and the initial interval is [0, 1].

The following are the codes:

Python

from __future__ import division, print_function
from numpy import log2, ceil, abs, cos

def bisection(fun, a, b, xtol=1e-6, ftol=1e-12):
    if fun(a) * fun(b) > 0:
        c = None
        msg = "The function should have a sign change in the interval."
    else:
        nmax = int(ceil(log2((b - a)/xtol)))
        for cont in range(nmax):
            c = 0.5*(a + b)
            if abs(fun(c)) < ftol:
                msg = "Root found with desired accuracy."
                break
            elif fun(a) * fun(c) < 0:
                b = c
            elif fun(b) * fun(c) < 0:
                a = c
            msg = "Maximum number of iterations reached."
    return c, msg

def fun(x):
    return cos(x) - x**2

print(bisection(fun, 0.0, 1.0))

With result

(0.824131965637207, 'Maximum number of iterations reached.')

Julia

function bisection(fun, a, b, xtol=1e-6, ftol=1e-12)
    if fun(a) * fun(b) > 0
        c = nothing
        msg = "The function should have a sign change in the interval."
    else
        nmax = ceil(log2((b - a)/xtol))
        for cont = 1:nmax
            c = 0.5*(a + b)
            if abs(fun(c)) < ftol
                msg = "Root found with desired accuracy."
                break
            elseif fun(a) * fun(c) < 0
                b = c
            elseif fun(b) * fun(c) < 0
                a = c
            end
            msg = "Maximum number of iterations reached."
        end
    end
    return c, msg
end

function fun(x)
    return cos(x) - x^2
end

println(bisection(fun, 0.0, 1.0))

With result

(0.824131965637207, "Maximum number of iterations reached.")

In this case, both codes are really close to each other. The Python code has 25 lines, while the Julia one has 27. As expected, the results given by them are the same.

Edit (2017-10-02)

As suggested by Edward Villegas, I decided to compare execution times. I used %timeit for IPython and @benchmark (from BenchmarkTools) for Julia.

In IPython, we have

%timeit bisection(fun, 0.0, 2.0)

with result

10000 loops, best of 3: 107 µs per loop

And in Julia, we have

@benchmark bisection(fun, 0.0, 2.0)

with result

BenchmarkTools.Trial:
  memory estimate:  48 bytes
  allocs estimate:  2
  --------------
  minimum time:     1.505 μs (0.00% GC)
  median time:      1.548 μs (0.00% GC)
  mean time:        1.604 μs (0.00% GC)
  maximum time:     38.425 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     10

And it seems that the Julia version is around 100 times faster than the Python counterpart.

Solution of the Schrödinger equation using Ritz method

In this post, we describe the solution of the Schrödinger equation using Ritz method and Hermite functions basis. This basis seems to be a good choice for the 1D Schrödinger equation since its an orthogonal basis over \((-\infty, \infty)\).

Transforming the equation to an algebraic one

We want so solve the time-independent Schrödinger equation

\begin{equation*} \left[-\frac{1}{2}\nabla^2 + V(x)\right] \psi = E\psi\, , \end{equation*}

where we are using natural units since they are a good choice for numeric calculations.

Solving equation is equivalent to solve the following variational equation

\begin{equation*} \delta \Pi[\psi] = 0\, , \end{equation*}

with

\begin{equation*} \Pi[\psi] = \frac{1}{2} \langle \nabla \psi, \nabla\psi\rangle + \langle \psi, V(x) \psi\rangle - E\langle \psi, \psi\rangle \, , \end{equation*}

being \(\psi\) the wave function, \(V(x)\) the potential and \(E\) the energy. This variational formulation is equivalent to the time-independent Schrödinger equation, and \(E\) works as a Lagrange multiplier to enforce that the probability over the whole domain is 1.

We can expand the wave function in an orthonormal basis, namely

\begin{equation*} \psi = \sum_{n=0}^N c_n u_n(x)\, , \end{equation*}

where \(u_n(x) \equiv \mu_n H_n(x) e^{-x^2/2}\) is a normalized Hermite function, \(\mu_n\) is the inverse of magnitude of the \(n\)th Hermite polynomial

\begin{equation*} \mu_n = \frac{1}{\sqrt{\pi^{1/2} n! 2^n}}\, , \end{equation*}

and \(c_n\) are the coefficients of the expansion. This representation is exact in the limit \(N \rightarrow \infty\).

If we replace the expansion in the functional, we obtain

\begin{equation*} \Pi_N = \sum_{m=0}^N\sum_{n=0}^N c_m c_n\left[ \frac{1}{2} \langle \nabla u_m, \nabla u_n\rangle + \langle u_m, V(x) u_n\rangle - E^N \delta_{mn}\right]\, . \end{equation*}

The integrand of the integral involving the two derivatives reads

\begin{align*} u_m' u_n' =& \left[2m \frac{\mu_{m-1}}{\mu_m}u_{m-1} - x u_m\right] \left[2n \frac{\mu_{n-1}}{\mu_n}u_{n-1} - x u_n\right]\\ =& 4mn\frac{\mu_{m-1} \mu_{n-1}}{\mu_m \mu_n} u_{n-1} u_{m-1} - 2m\frac{\mu_{m-1}}{\mu_{m}}x u_{m-1} u_n\\ &- 2n\frac{\mu_{n-1}}{\mu_{n}}x u_{n-1} u_m + x^2 u_m u_n \end{align*}

Thus, the kinetic energy term reads

\begin{align*} \langle \nabla u_m, \nabla u_n \rangle =& 4mn\frac{\mu_{m-1} \mu_{n-1}}{\mu_m \mu_n} \langle u_{n-1}, u_{m-1}\rangle - 2m\frac{\mu_{m-1}}{\mu_{m}} \langle u_{m-1}, x u_n\rangle\\ &- 2n\frac{\mu_{n-1}}{\mu_{n}} \langle u_{m}, x u_{n - 1}\rangle + \langle u_m, x^2 u_n\rangle\\ =& 4mn \frac{\mu_{m-1}^2}{\mu_m^2}\delta_{mn} - 2m \frac{\mu_{m-1}}{\mu_m} \alpha_{m-1, n} - 2n \frac{\mu_{n-1}}{\mu_n} \alpha_{m, n-1} + \beta_{mn} \, , \end{align*}

with

\begin{equation*} \alpha_{m,n} \equiv \langle u_{m}, x u_n\rangle = \begin{cases} \sqrt{\frac{n + 1}{2}} & m=n +1\\ \sqrt{\frac{n}{2}} & m=n - 1\\ 0 & \text{otherwise}\end{cases}\, , \end{equation*}

and

\begin{equation*} \beta_{m,n} \equiv \langle u_{m}, x^2 u_n\rangle = \begin{cases} \frac{\sqrt{n(n-1)}}{2} & m = n - 2\\ \frac{2n + 1}{2} & m = n \\ \frac{\sqrt{(n+1)(n+1)}}{2} & m = n + 2 \\ 0 & \text{otherwise}\end{cases}\, . \end{equation*}

The functional is rewritten as

\begin{align*} \Pi_N =& \sum_{m=0}^N \sum_{n=0}^N c_m c_n \left[2mn \frac{\mu^2_{m-1}}{\mu^2_m}\delta_{mn} - m\frac{\mu_{m-1}}{\mu_m}\alpha_{m - 1, n} - n\frac{\mu_{n-1}}{\mu_n}\delta_{m, n-1} \right.\nonumber \\ &\left. - \frac{1}{2}\beta_{mn} + \langle u_m, V(x) u_n\rangle - E^N \delta_{mn}\right] \, . \end{align*}

Taking the variation

\begin{equation*} \delta \Pi_N = 0\, , \end{equation*}

that in this case is equivalent to

\begin{equation*} \frac{\partial \Pi_N}{\partial c_i} = 0\quad \forall i=0, 1, \cdots N\, , \end{equation*}

yields to

\begin{equation*} [H]\lbrace c\rbrace = E^N\lbrace c\rbrace \, , \end{equation*}

where the components of the matrix \([H]\) are given by

\begin{equation*} H_{mn} = 2mn \frac{\mu^2_{m-1}}{\mu^2_m}\delta_{mn} - m\frac{\mu_{m-1}}{\mu_m}\alpha_{m - 1, n} - n\frac{\mu_{n-1}}{\mu_n}\delta_{m, n-1} - \frac{1}{2}\beta_{mn} + \langle u_m, V(x) u_n\rangle\, . \end{equation*}

The last integral can be computed using Gauss-Hermite quadrature. And we will need more Gauss points if we want to integrate higher-order polynomials. This method would work fine for functions that can be approximated by polynomials.

Examples

A Python implementation of this method is presented in this repo.

For all the examples we use the following imports

from __future__ import division, print_function
import numpy as np
from hermite_QM import *

Quantum harmonic oscilator

In this case the (normalized) potential is given by

\begin{equation*} V(x) = \frac{1}{2} x^2 \end{equation*}

and the exact normalized eigenvalues are given by

\begin{equation*} E_n = n + \frac{1}{2} \end{equation*}

The following snippet computes the first 10 eigenvalues and plot the corresponding eigenstates

potential = lambda x: 0.5*x**2
vals, vecs = eigenstates(potential, nterms=11, ngpts=12)
print(np.round(vals[:10], 6))
fig, ax = plt.subplots(1, 1)
plot_eigenstates(vals, vecs, potential);

with results

[ 0.5  1.5  2.5  3.5  4.5  5.5  6.5  7.5  8.5  9.5]
/images/hermite_ritz_harmonic.svg

Absolute value potential

potential = lambda x: np.abs(x)
vals, vecs = eigenstates(potential)
print(np.round(vals[:10], 6))
fig, ax = plt.subplots(1, 1)
plot_eigenstates(vals, vecs, potential, lims=(-8, 8));

with results

[ 0.811203  1.855725  2.57894   3.244576  3.826353  4.38189   4.895365
  5.396614  5.911591  6.421015]
/images/hermite_ritz_abs.svg

Cubic potential

potential = lambda x: np.abs(x/3)**3
vals, vecs = eigenstates(potential)
print(np.round(vals[:10], 6))
fig, ax = plt.subplots(1, 1)
plot_eigenstates(vals, vecs, potential, lims=(-6, 6));

with results

[ 0.180588  0.609153  1.124594  1.681002  2.272087  2.889805  3.530901
  4.191962  4.871133  5.566413]
/images/hermite_ritz_cubic.svg

Harmonic with quartic perturbation

potential = lambda x: 0.5*x**2 + 0.1*x**4
vals, vecs = eigenstates(potential, nterms=20, ngpts=22)
print(np.round(vals[:10], 6))
fig, ax = plt.subplots(1, 1)
plot_eigenstates(vals, vecs, potential, lims=(-5, 5))

with results

[  0.559146   1.769503   3.138624   4.628884   6.220303   7.899789
   9.658703  11.489094  13.38638   15.361055]
/images/hermite_ritz_pert_harm.svg

A Jupyter Notebook with the examples can be found here.

Using style sheets with matplotlib

We can format plots easily using the style package in matplotlib all. The main idea is to create a file with some of the parameters that want to be defined (that can also be accessed through rcParams).

This post is not a tutorial on how to use those, for that you can check the style sheet reference. Here, I just want to play with some of these parameters to create three different styles. The first two examples present the style of an (infamous by some) software, that is probably used for most people as their visualization platform, while the third one is just a clean style. All the files used here can be download here.

For all the examples below the following imports are done:

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

First example: MS 2003

In our first example we want to reproduce the style that we used to see more than a decade ago as default.

The following is the content of the file MS2003.mplstyle

font.family : sans-serif

axes.facecolor : c0c0c0
axes.edgecolor : black
axes.prop_cycle : cycler('color',['000080', 'FF00FF', 'FFFF00', '00FFFF','800080', '800000', '008080', '0000FF'])
axes.grid : True

axes.spines.left   : True
axes.spines.bottom : True
axes.spines.top    : True
axes.spines.right  : True

grid.color : black
grid.linestyle : -

lines.linewidth : 1

figure.figsize : 5, 3

legend.fancybox : False
legend.frameon : True
legend.facecolor : white
legend.edgecolor : black
legend.loc : center left

The following code use this style

style = "MS2003.mplstyle"
with plt.style.context(style):
    x = np.linspace(0, 4, 100)
    y = np.sin(np.pi*x + 1e-6)/(np.pi*x + 1e-6)
    fig = plt.figure()
    ax = plt.subplot(111)
    for cont in range(5):
        plt.plot(x, y/(cont + 1), label=cont)

    plt.gca().xaxis.grid(False)
    box = ax.get_position()
    ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
    plt.legend(bbox_to_anchor=(1, 0.5))

and this is the result

/images/MS2003_style.svg

Second example: MS 2007

In the second example we want to reproduce the offspring of the style in the first example. This is definitely an improvement over the previous style, but it is not perfect.

The following is the content of the file MS2007.mplstyle

font.family : sans-serif

axes.facecolor : white
axes.edgecolor : 4d4d4d
axes.prop_cycle : cycler('color',['4573a7', 'aa4644', '89a54e', '71588f','4298af', 'db843d', '93a9d0', 'd09392'])
axes.grid : True
axes.linewidth : 0.5

axes.spines.left   : True
axes.spines.bottom : True
axes.spines.top    : False
axes.spines.right  : False

lines.linewidth : 2

grid.color : 4d4d4d
grid.linestyle : -
grid.linewidth : 0.5

figure.figsize : 5, 3

legend.fancybox : False
legend.frameon : False
legend.facecolor : white
legend.edgecolor : 4d4d4d
legend.loc : center left

The following code use this style

style = "MS2007.mplstyle"
with plt.style.context(style):
    x = np.linspace(0, 4, 100)
    y = np.sin(np.pi*x + 1e-6)/(np.pi*x + 1e-6)
    fig = plt.figure()
    ax = plt.subplot(111)
    for cont in range(5):
        plt.plot(x, y/(cont + 1), label=cont)

    plt.gca().xaxis.grid(False)
    box = ax.get_position()
    ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
    plt.legend(bbox_to_anchor=(1, 0.5))

and this is the result

/images/MS2007_style.svg

Third example: a clean style

The last example is a clean style that uses a color palette taken from ColorBrewer.

The following is the content of the file clean_style.mplstyle

font.family : sans-serif

axes.facecolor : white
axes.prop_cycle : cycler('color',['e41a1c', '377eb8', '4daf4a', '984ea3', 'ff7f00', 'ffff33', 'a65628', 'f781bf'])
axes.linewidth : 0.0
axes.grid : True

lines.linewidth : 1.5

xtick.direction : in
ytick.direction : in

grid.color : c7dedf
grid.linestyle : -
grid.linewidth : 0.3

figure.figsize : 6, 4

legend.fancybox : False
legend.frameon : False
legend.loc : best

The following code use this style

style = "clean.mplstyle"
with plt.style.context(style):
    x = np.linspace(0, 4, 100)
    y = np.sin(np.pi*x + 1e-6)/(np.pi*x + 1e-6)
    fig = plt.figure()
    ax = plt.subplot(111)
    for cont in range(5):
        plt.plot(x, y/(cont + 1), label=cont)

    plt.legend()

and this is the result

/images/clean_style.svg

We can also use files that are stored remotely. For example, in the third example we could have used the following URL:

style = "https://raw.githubusercontent.com/nicoguaro/matplotlib_styles/master/styles/clean.mplstyle"

Resources

As I mentioned above, the objective of this post was jut to create some simple-enough style-sheets for matplotlib and see them in action.

This post does not pretend to be a guide for good plots/visualization. For that matter you better look for some references, for example:

Also, I found really useful the following tools:

  • ColorBrewer2 allows to pick colormaps with different criteria (quantitative/qualitative, printer friendly, colorblind friendly).
  • ColRD has plenty of color palettes. It also has the option to generate palettes from images.
  • Colorgorical is a tool to make categorical color palettes for information visualization.

You can find the snippets present in this post in this Jupyter notebook.

Piecewise Hermite interpolation

In this post, we find the Hermite interpolation functions for the domain [-1, 1]. And then, we use it for a pieciwise interpolation. Notice that this interpolation has \(C^1\) continuity compared to the \(C^0\) continuity that is common in Lagrange interpolation.

To compute the polynomials explicitly we use sympy.

from __future__ import division
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

We want to find a set of basis function that satisfy the following

\begin{align*} N_1(x_1) &= 1\\ N_1(x_2) &= 0\\ N_2(x_1) &= 1\\ N_2(x_2) &= 0\\ N'_3(x_1) &= 1\\ N'_3(x_2) &= 0\\ N'_4(x_1) &= 1\\ N'_4(x_2) &= 0 \end{align*}

These can be written as

\begin{equation*} \begin{bmatrix} 1 &x_1 &x_1^2 &x_1^3\\ 1 &x_2 &x_2^2 &x_2^3\\ 0 &1 &2x_1 &3x_1^2\\ 0 &1 &2x_2 &3x_2^2 \end{bmatrix} \begin{bmatrix} a_{11} &a_{12} &a_{13} &a_{14}\\ a_{21} &a_{22} &a_{23} &a_{24}\\ a_{31} &a_{32} &a_{33} &a_{34}\\ a_{41} &a_{42} &a_{43} &a_{44} \end{bmatrix} = \begin{bmatrix} 1 &0 &0 &0\\ 0 &1 &0 &0\\ 0 &0 &1 &0\\ 0 &0 &0 &1 \end{bmatrix} \end{equation*}

The formula for interpolation would be

\begin{equation*} f(x) \approx N_1(x) u_1 + N_2(x) u_2 + |J|(N_3(x) u'_3 + N_4(x) u'_4)\quad \forall x\in [a, b] \end{equation*}

with \(|J| = (b - a)/2\) the Jacobian determinant of the transformation.

x1, x2, x = sym.symbols("x1 x2 x")
V = sym.Matrix([
    [1, x1, x1**2, x1**3],
    [1, x2, x2**2, x2**3],
    [0, 1, 2*x1, 3*x1**2],
    [0, 1, 2*x2, 3*x2**2]])
V
\begin{equation*} \left[\begin{matrix}1 & x_{1} & x_{1}^{2} & x_{1}^{3}\\ 1 & x_{2} & x_{2}^{2} & x_{2}^{3}\\ 0 & 1 & 2 x_{1} & 3 x_{1}^{2}\\ 0 & 1 & 2 x_{2} & 3 x_{2}^{2}\end{matrix}\right] \end{equation*}

We can see that the previous matrix is a confluent Vandermonder matrix [1]. It is similar to a Vandermonde matrix for the first \(n\) nodes and the derivatives of each row for the following ones.

The coefficients for our polynomials are given by the inverse of this matrix.

sym.simplify(V.inv())
\begin{equation*} \left[\begin{matrix}\frac{x_{2}^{2} \left(3 x_{1} - x_{2}\right)}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & \frac{x_{1}^{2} \left(x_{1} - 3 x_{2}\right)}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & - \frac{x_{1} x_{2}^{2}}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}} & - \frac{x_{1}^{2} x_{2}}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}}\\- \frac{6 x_{1} x_{2}}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & \frac{6 x_{1} x_{2}}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & \frac{x_{2} \left(2 x_{1} + x_{2}\right)}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}} & \frac{x_{1} \left(x_{1} + 2 x_{2}\right)}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}}\\\frac{3 x_{1} + 3 x_{2}}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & - \frac{3 x_{1} + 3 x_{2}}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & - \frac{x_{1} + 2 x_{2}}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}} & - \frac{2 x_{1} + x_{2}}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}}\\- \frac{2}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & \frac{2}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & \frac{1}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}} & \frac{1}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}}\end{matrix}\right] \end{equation*}

And we replace the values for the extreme nodes, -1 and 1.

V_inv = sym.simplify(V.subs({x1:-1, x2:1}).inv())
V_inv
\begin{equation*} \left[\begin{matrix}\frac{1}{2} & \frac{1}{2} & \frac{1}{4} & - \frac{1}{4}\\ - \frac{3}{4} & \frac{3}{4} & - \frac{1}{4} & - \frac{1}{4}\\ 0 & 0 & - \frac{1}{4} & \frac{1}{4}\\ \frac{1}{4} & - \frac{1}{4} & \frac{1}{4} & \frac{1}{4}\end{matrix}\right] \end{equation*}

The polynomials are the product of the coefficients and the monomials

sym.factor(V_inv.T * sym.Matrix([1, x, x**2, x**3]))
\begin{equation*} \left[\begin{matrix}\frac{1}{4} \left(x - 1\right)^{2} \left(x + 2\right)\\- \frac{1}{4} \left(x - 2\right) \left(x + 1\right)^{2}\\\frac{1}{4} \left(x - 1\right)^{2} \left(x + 1\right)\\\frac{1}{4} \left(x - 1\right) \left(x + 1\right)^{2}\end{matrix}\right] \end{equation*}

The interpolation basis would be

\begin{align*} N_1 (x) &= \frac{1}{4} (x - 1)^2 (2 + x)\\ N_2 (x) &= \frac{1}{4} (x + 1)^2 (2 - x)\\ N_3 (x) &= \frac{1}{4} (x - 1)^2 (x + 1)\\ N_4 (x) &= \frac{1}{4} (x + 1)^2 (x - 1)\, , \end{align*}

and the following function computes the interpolation for a given function and derivative

def hermite_interp(fun, grad, x0=-1, x1=1, npts=101):
    jaco = (x1 - x0)/2
    x = np.linspace(-1, 1, npts)
    f1 = fun(x0)
    f2 = fun(x1)
    g1 = grad(x0)
    g2 = grad(x1)
    N1 = 1/4*(x - 1)**2 * (2 + x)
    N2 = 1/4*(x + 1)**2 * (2 - x)
    N3 = 1/4*(x - 1)**2 * (x + 1)
    N4 = 1/4*(x + 1)**2 * (x - 1)
    interp = N1*f1 + N2*f2 + jaco*(N3*g1 + N4*g2)
    return interp

In this case, we interpolate the sinc function

def fun(x):
    return np.sin(2*np.pi*x)/(2*np.pi*x)


def grad(x):
    return np.cos(2*np.pi*x)/x - np.sin(2*np.pi*x)/(2*np.pi*x**2)

The following snippet compute the interpolation and plot it.

a = 2
b = 5
nels = 7
npts = 200
x = np.linspace(a, b, npts)
y = fun(x)
plt.plot(x, y, color="black")
xi = np.linspace(a, b, num=nels, endpoint=False)
dx = xi[1] - xi[0]
for x0 in xi:
    x1 = x0 + dx
    x = np.linspace(x0, x1, npts)
    y = hermite_interp(fun, grad, x0=x0, x1=x1, npts=npts)
plt.plot(x, y, linestyle="dashed", color="#4daf4a")
plt.plot([x[0], x[-1]], [y[0], y[-1]], marker="o", color="#4daf4a",
         linewidth=0)
plt.xlabel("x")
plt.ylabel("y")
plt.legend(["Exact function", "Interpolation"])
plt.show()
/images/sinc_interp.svg

References

[1] Walter Gautschi (1962). On inverses of Vandermonde and confluent Vandermonde matrices. Numerische Mathematik, 4 117-123.

Herramientas (auxiliares) para Investigación

En esta entrada quiero hablar de algunas herramientas (...o algo así) que son útiles en el día a día de la investigación, pero que no suelen ser tan populares por ser en cierta forma tangenciales al área específica en la que cada uno trabaja.

Scripting

Un script  es un programa usualmente simple que realiza una serie de tareas. ¡Lo siento! trato de ceñirme al castellano pero no me gustan las traducciones para este término (archivo de órdenes, archivo de procesamiento por lotes o guion). Si existe una tarea que debe ser hecha más de... cinco veces (el número varía de acuerdo a la paciencia), entonces creo que es algo que podemos pedirle al computador que haga esta tarea por nosotros. En otras palabras: podemos automatizar ese trabajo. Algunas tareas que pueden ser una buena alternativa para automatizar son: renombrar 100 archivos, convertir un archivo de un formato a otro (p. ej., STL a OBJ), leer 387 archivos con información sobre el clima y graficar su evolución temporal (temperatura mínima, máxima y promedio). Estas tareas pueden ser fáciles de hacer a mano, pero por la cantidad de trabajo que involucran son tediosas.

Lo primero es hacernos con un lenguaje de scripting. Algunas opciones son Python, Bash, Julia, Matlab/Octave, Scilab. Dejándome llevar por el sesgo, recomiendo usar Python.

Gráficos y esquemas

Una imagen vale más que mil palabras, o eso dice el dicho. Personalmente, me parece absolutamente cierto y trato de hacer dibujitos para poder entender mejor algo o lograrlo explicar mejor. Lo primero que me gustaría mencionar es la diferencia entre imágenes de mapas de bits (o ráster) e imágenes vectoriales.

  • Imagen de mapa de bits: es una imagen que está representada por un arreglo (o rejilla rectangular) de pixeles. Dicho de otro modo, se almacena la información de color que hay en cada punto de la imagen. Los formatos más populares almacenan la información comprimida. Para gráficos con alto contraste (como esquemas o diagramas) el mejor formato es PNG. Si se tiene una animación, GIF sería preferible. Y para el caso de fotografías es mejor utilizar JPG.
  • Imagen vectorial: es una imagen que está formada por entidades geométricas. En esta no se almacena la información punto a punto sino la construcción de las formas que la constituyen. Por esta razón, estas imágenes no se pixelan pues la información que se tiene es de cómo construirla. Este tipo de imágenes es la mejor opciones para esquemas y diagramas, pues sólo se almacena la información de los trazos y texto que se añadan en ellas (ver Figura 1). El estándar de facto para este tipo de imágenes es PDF —es el que suelo para incluir en mis documentos \(\LaTeX\), aunque existe forma de embeber SVG en \(\LaTeX\) pero es algo que aún no exploro—. A pesar de que PDF sea el estándar, el formato a preferir es SVG (Scalable Vector Graphics) que es un estándar a través de la internet y la mayoría de navegadores modernos permiten visualizar.

Recapitulando, deberíamos usar imágenes JPG para fotografías y SVG para esquemas/diagramas. Otro atributo que puede ser de utilidad es el manejo de capas, SVG permite esto.... y de formatos ráster tenemos la opción de usar TIFF.

En cuanto a software para generar/editar este tipo de imágenes debo decir que existen gran cantidad de programas que permiten exportar a estos formatos: Python/Matplotlib, Matlab, Inkscape, Adobe Illustrator, GIMP, Photoshop, LibreOffice. Si el gráfico es generado a partir de un cálculo o de una serie de datos yo uso Matplotlib. Si en cambio, queremos hacer un esquema o —como suelo decir— un muñeco mi herramienta favorita es Inkscape. Este programa pretende ser una alternativa libre a programas como Adobe Illustrator —y sí que lo consigue. Obviamente, se podría usar Illustrator o Corel Draw para esta tarea. Si el único uso sería hacer esquemas técnico, creo que sería un desperdicio.

/images/sensor_hall.png
Figura 1. Esquema de sensado de campo magnético de imanes permanentes por efecto Hall. Sacado de mi trabajo de grado de Ingeniería Física.

Tomar notas

Supongo que para algunos parecería un poco trivial hablar de "tomar notas" y mucho más viniendo de alguien que no tenía cuadernos en el bachillerato, pero como soy un poco terco creo que igual escribiré un poco sobre esto. Lo primero que me gustaría mencionar es que recuerdo que me hablaran de este tema en el colegio, pero nunca hubo nada formal respecto a desarrollar estas habilidades. Buscando en la web, hay gran cantidad de información. Incluso el artículo de Wikipedia en inglés es interesante. No hay nada mejor para escribir que tener una pluma y un papel de buen gramaje, por eso es que aún uso un cuaderno en donde llevo cuenta de lo que hago en mi investigación y tomo notas. Sin embargo, este esquema es bastante lineal y deja por fuera oportunidades más contemporáneas. Es decir ¿por qué conformarse con un documento en este tiempo de hiper-documentos? Las ventajas de tomar notas digitalmente saltan a la vista, en un hiper-documento se pueden tener enlaces, embeber imágenes, video y sonido.

En cuanto a herramientas, acá incluyo una pequeña lista:

  • Evernote es probablemente la herramienta más popular para tomar notas. Es multiplataforma, Freemium (funcionalidad básica gratuita y avanzada paga), y cuenta con muchas opciones. Yo la uso, pero no mucho en mi investigación.
  • Zim es una wiki offline. Tiene gran cantidad de opciones como calendario, ecuaciones con código \(\LaTeX\), imágenes... en fin. El pero que le encuentro es que no he logrado configurar las ecuaciones en Windows (y en mi oficina debo usar Ruindows :-/).
  • Docear esta es una herramienta pensada, principalmente, para manejar bibliografía. Sin embargo, permite tomar notas y, en general, manejar la información de la investigación. La característica más (o menos) atractiva es que funciona alrededor de mapas mentales.
  • Zotero también es una herramienta para manejar bibliografía, aunque permite manejar algo de toma de notas (al menos alrededor de la bibliografía).
  • Mendeley es muy parecido al anterior, aunque con más funcionalidades. El pero más grande que le encuentro es que en 2013 fue comprado por Elsevier.

En cuanto a manejo de bibliografía también me gustaría mencionar EndNote que es el programa con mayor trayectoria y  JabRefque es el que he usado por más tiempo. Algunas referencias interesantes comparando manejadores de bibliografía están acá: [A] [B] [C].

Reconstrucción de gráficos

Es común encontrarse con información presentada en forma de gráficos. También es común que queramos tener los datos numéricos de estos gráficos para poder compararlos con los nuestros. Para saber si nuestras medidas/simulaciones/métodos dan resultados similares a otros presentados en la literatura. Para esta tarea se pueden usar potentes software de procesamiento de imágenes, u otros más modestos diseñados específicamente para este fin.

/images/pointplot.jpg
Figura 2. Gráfico original.
/images/digitized_pts.png
Figura 3. Gráfico procesado en Engauge Digitizer. Algunos puntos fueron seleccionados (automáticamente) para obtener sus coordenadas.
  • Digitizer of XY chart este es un plugin para Libreoffice/OpenOffice y exporta el resultado a la hoja de cálculo actual, es simple y fácil de usar.
  • Engauge Digitizer, es el que normalmente uso cuando necesito hacer esta tarea (ver Figuras XX). Es de código abierto (y libre) y tiene una buena cantidad de opciones para facilitar la tarea.
  • Plot Digitizer no tengo mucha información sobre este (pues nunca lo he usado), salvo que está escrito en Java.
  • ImageJ este es un (completo) programa de procesamiento de imágenes que está escrito en Java. No lo he usado para esta tarea de manera regular, pero podría usarse para ello.

Visualización científica

La visualización científica que se encarga de generar gráficos que permitan visualizar "datos científicos" para facilitar el entendimiento que hay detrás de los datos. Para esta labor muchos hemos usado lenguajes de scripting como Matlab/Octave, Scilab o Python (con Matplotlib o Mayavi). Sin embargo, como la visualización se trata de algo visual —¿como si no?—, es bueno contar con una herramienta que permita generar y cambiar los gráficos de manera interactiva, aunque siempre debemos automatizar la mayor cantidad de trabajo posible (la pereza siempre ha sido uno de los móviles más grandes de la humanidad, hay que aceptarlo).

  • MayaVi, este es un programa escrito en Python que usa VTK. Es una herramienta muy versátil y la gran ventaja que tiene es que puede usarse dentro de scripts de Python.
  • Paraview, este programa también está basado en VTK y permite paralelizar las labores (para los computadores con múltiple núcleo y los clusters). Abajo incluyo un video generado en Paraview para mostrar sus capacidades.
  • Visit, este programa también está basado en VTK, nunca lo he usado pero quise incluirlo porque la gente dice que puede ser más intuitivo que Paraview.
  • Tecplot, este programa es muy popular en Purdue. Creo que inicialmente fue pensado para CFD, pero se ha expandido mucho. En cuanto a gráficos 3D no me parece mejor que ParaView, sin embargo, las capacidades de gráficos 2D (gráficos XY, y demás) lo hacen atractivo.
  • Scavis, este está escrito en Java. No lo conocía hasta que inicie la escritura de esta entrada pero me llamó la atención y quise incluirlo en la lista. Algo que me llagmó la atención es que permite scripting en varios lenguajes: Java, Python, Ruby, BeanShell y Matlab/Octave.
  • Origin, nunca lo he usado pero no lo quise dejar por fuera porque siempre he oído hablar maravillas de él (probablemente compárandolo con Excel... pero no puedo opinar).

Simulación de fuego cruzado en Vimeo.

Manejo de versiones

El manejo de versiones es la administración de cambios en documentos, código fuente y otro tipo de información. Esto puede hacerse de forma manual, pero es fácil cometer errores o remplazar la versión de un código fácilmente, y por esto es recomendable usar un software que facilite el trabajo. La idea es tener un lugar (repositorio) en donde se almacenan las versiones y los cambios, y llevar un registro de estos. De esta forma se puede volver a una versión anterior de los documentos y varias personas pueden trabajar conjuntamente.  Existen dos paradigmas (o arquitecturas) para el manejo de versiones: centralizada y distribuida. En la primera existe un repositorio centralizado en donde se encuentra toda la información. En la arquitectura distribuida cada usuario tiene una copia del respositorio. Personalmente sólo he usado Git, que entra en la categoría distribuida y es uno de los software de manejo de versiones más populares actualmente; lo usan compañías como Google, Facebook y Netflix.

Un ejemplo puede verse en este repositorio, en donde está el documento de trabajo de grado de Santiago Echeverri, el cual tuve la oportunidad de asesorar. Este documento lo editamos conjuntamente mientras él estaba en Medellín  y yo me encontraba en Estados Unidos. El documento se hizo en el lenguaje de marcadores \(\LaTeX\).

Además de tener un control sobre las versiones y poder acceder a versiones anteriores, es útil poder almacenar la información en un lugar asequible desde cualquier lugar del mundo con una conexión a internet. Esto puede lograrse con un servidor propio, obviamente, o también a través de un proveedor externo. Dos proyectos que  son muy populares para alojar repositorios y manejar sus versiones son (comparación entre Github y BitBucket):

  • Github  es el más popular en este momento. Permite tener proyectos con un número ilimitado de colaboradores. Para tener un repositorio privado es necesario pagar.
  • BitBucket la principal ventaja es que permite tener repositorios privados sin la necesidad de pagar. Sólo es gratuito para proyectos con 5 colaboradores o menos (o para proyectos académicos).

Enlaces sugeridos

  1. Software Carpentry. http://software-carpentry.org/
  2. Python Scientific Lecture Notes. https://scipy-lectures.github.io/

Seguro dejé mucho temas por fuera así como herramientas dentro de algún tópico. Si ese es el caso, agradecería que me lo digan en los comentarios.

La zorra y las uvas

Mirando en internet no encontré una versión en español de la fábula de Esopo "La zorra y las uvas" que me gustara (o que al menos rimara :P). De hecho, el artículo de Wikipedia en español es bastante más malo que su equivalente en inglés.  Por esto, y por petición de una amiga, hice mi versión con rima asonante. Y reza así

La Zorra y las uvas

La zorra que con ansias uvas buscaba
no pudo alcanzarlas por lo alto que estaban.
Giró y exclamó, con risa lastimera:
¡Están verdes! No valen mi espera.

Tomando notas de clase con mi tableta

El año pasado el niño dios me regaló una Samsung Galaxy Note 10.1 (modelo 2012). Esta tableta viene diseñada para tomar notas, pues trae un S-Pen (fabricado por Wacom, que creo que fabrican la tecnología para la mayoría de tabletas comerciales) y varias aplicaciones pensadas para ello. Como: S Note, PhotoShop Touch, Crayon Physics (tienen que jugarlo es muy divertido).

Decidí entonces hacer un experimento tomando notas de clase con mi tableta completamente en las dos materias que cursé en el "semestre" de primavera de 2014. En Optimización usé el lapiz que viene con la tableta y para Física Molecular usé un Bamboo Stylus que, se supone, es diseñado específicamente para esta tableta. Las notas las tomé con el software por defecto (S Note) y convertí a PDF (las notas fueron hechas enteramente con el software):

Aunque la sensación al tacto del Bamboo Stylus es mucho más agradable, el desempeño del S Pen es mucho mejor. En ambos casos, la sensación es muy diferente a la de usar papel y lápiz (lapicero/pluma), en gran medida porque la fricción es muy poca en comparación (algo que es mejor con los lapices conductores). La experiencia de usuario del lápiz varía mucho entre una aplicación y otra, lo que genera buenas expectativas. En cuanto a tomar notas: es muy recomendado para tomar notas cortas (como reuniones); no tanto para tomar notas largas (como notas de clase), a menos que se quiera tener una copia digital fácilmente.

Apps

La tecnología detrás de la tableta es la parte más importante. De manera simple y rápida, el lapiz tiene una antena y la tableta una rejilla de antenas que están sintonizadas en la misma frecuencia (531 kHz, en este enlace explican su funcionamiento un poco mejor). Sin embargo, también es importante tener en cuenta el software que se usa para tomar nota:

S Note

  • Es la aplicación nativa para tomar notas.
  • Tiene reconocimiento de formas y ecuaciones. Las ecuaciones deben ser simples, no esperen tomar notas de Mecánica Cuántica o Mecánica del Medio Continuo usado esta herramienta.
  • Permite insertar imágenes y grabar sonidos dentro de la aplicación.
  • Tiene "palm rejection".
  • El tamaño de papel no es configurable.

Papyrus

  • Tiene "palm rejection" y permite configurar el dedo como borrador (que puede ser útil).
  • Permite tener varios tamaños de papel y diferentes fondos.
  • Permite hacer rectángulos y elipses.
  • Los pinceles son muy básicos.

Note Anytime

  • Permite tener varios tamaños de papel.
  • Almacena los trazos vectorialmente y se pueden editar (cambiar color, tamaño o estilo posteriormente).
  • Tiene varios pinceles y además son configurables.
  • No puede convertir las notas en PDF.
  • No tiene "palm rejection".

Quill

  • Es open source, aunque si lo adquieren a través de la "Google Store" tiene costo.
  • El algoritmo de reconocimiento de trazos es mejor que el de los otros software (y es ajustable).
  • Aún está en una etapa preliminar.

Lástimosamente no hay una versión de Paper (de FiftyThree) para Android, y, aunque existe una versión de Bamboo Paper (de Wacom) para Android no es compatible con esta tableta.

En conclusión

Remplazar la toma de notas "tradicionales" por notas  "digitales" puede ser un poco apresurado, pero es realizable. La aplicación que viene por defecto (S Note) permite realizar el trabajo "out of the box", pero podría ser mejor. Las aplicaciones que me parecen más promisorias son: Note Anytime y Quill.

Palabras terminadas en "ar", "er" o "ir"

Recuerdo cuando estaba en el colegio que en ocasiones mencionaban las palabras terminadas en "ar", "er" o "ir" —iguales al final de nuestros verbos en infinitivo—. Eran pocos ejemplos los que mencionaban, así que algunos años después me dio el arrebato de hacer una lista… y hela aquí:

AR AR AR ER IR
agar espectacular paramilitar alfiler cachemir
ajuar estelar particular alquiler casimir
alar extracurricular peculiar amater elixir
albañar familiar pendular anteayer emir
albar fular peninsular antier faquir
alfar globular perpendicular atelier kefir
alizar hangar pilar ayer menhir
almiar hogar pinar bachiller nadir
almimbar ijar planar bereber porvenir
alminar impar polar brigadier sir
altar insular policelular canciller souvenir
alveolar intercelular pulgar cualquier visir
antealtar intraocular pulmonar doquier zafir
antisolar jaguar quasar dossier  
ar juglar radar enser  
auricular lagar radicular ester  
avatar lanar rectangular menester  
axilar lar reticular mercader  
azahar limonar samovar neceser  
azar llar secular palier  
bacillar lobular seglar placer  
bar lodazar semicircular plumier  
bazar lugar semilunar premier  
bicelular lunar similar primer  
bienestar malar solar rosicler  
bifilar maleolar subcelular salmer  
biliar malvar supraclavicular somier  
billar mandibular supralunar sumiller  
binocular manglar telar taller  
bipolar manillar testicular tercer  
broncopulmonar manjar tisular veguer  
bulevar mar trifilar    
calamar maxilar trufar    
calcañar melocotonar tubular    
capilar melonar tular    
capsular membrillar ultramar    
caviar miliar uvular    
celular militar vascular    
centenar millar vehicular    
ciliar molar vermicular    
clavicular mollar vesicular    
collar monocelular vulgar    
consular mular yugular    
coplanar muscular zar    
corpuscular nuclear      
cuadrangular ocular      
cular olivar      
curricular pajar      
dactilar paladar      
dinar palatar      
dispar palomar      
ejemplar papilar      
escolar par      

Respecto a la teletransportación en Star Trek

Alguna vez una amiga me preguntó sobre cuál era más fácil de dos métodos para teletransportar personas. Uno consistía en "desconfigurar" y "reconfigurar" a la persona (como información, similar a Start Trek y la otra en "crear un atajo" entre dos lugares del espacio-tiempo para pasar a la persona. Obviamente, en ambos casos suponiendo que tales cosas se pueden hacer.

Respecto a la creación del atajo en el espacio-tiempo, creo que no he estudiado las ecuaciones de Relatividad General para poder hacer las cuentas, aunque tal vez en Física Pasión nos regalen estos cálculos.

La primera de las alternativas consiste, básicamente, en la representación de una persona como información y su procesamiento de "desconfiguración", transmisión y procesamiento de "reconfiguración". Ya que los humanos estamos hechos, mayormente, de agua asumiremos una persona de agua para los cálculos. Además tomaremos como masa de referencia \(70\ \mbox{kg}\).

Algunos datos que usaremos son:

  • La masa molar del agua es \(M_{H_2O}=18,01528\ \mbox{g/mol}\);
  • El número de Avogadro es \(N_A=6,022 \times 10^{23}/\mbox{mol}\).

Entonces tendríamos que el número de moléculas en una persona es

\begin{equation*} m_\text{persona} = \frac{70\ \mbox{kg}\ \mbox{moles}} {18,015\times 10^{-3} \text{kg}} =3885,7\ \mbox{moles}\, , \end{equation*}

y el número de moléculas sería

\begin{equation*} N_\text{moleculas} = N_A m_\text{persona} = 2,34\times 10^{27}\ \mbox{moleculas}\, . \end{equation*}

El agua tiene 12 modos vibracionales y 6 traslacionales, y además 40 números cuánticos electrónicos. Lo que suma un total de 58 grados de libertad por cada molécula. Así que el número total de grados de libertad es

\begin{equation*} N_{GDL} = 1,36\times 10^{29}\, . \end{equation*}

En 2011 hubo un record de transmisión de 186 Gb/s, considerando esta tasa el tiempo que tardaría en transmitirse la totalidad de los datos (asumiendo datos representados como reales de 32 bits) tardaría \(2,3\times 10^{16}\ \mbox{s}\) o, equivalentemente, \(1,70\times 10^{10}\) años (la edad estimada del universo es de \(1,37\times 10^{10}\) años).

Si quisiéramos que esto pasara ealmente rápido, digamos en 1 milisegundo, la tasa de transmisión debería ser \(2,3\times 10^{19}\) más rápida que el récord que alcanzaron en el CERN… algo que es inconcebible para nosotros actualmente.