Skip to main content

Nicolás' blog

Numerical methods challenge: Day 22

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.

Finite Difference: Eigenproblems

Today we have a Finite difference method to compute the vibration modes of a beam. The equation of interest is

\begin{equation*} \nabla^ 4 u = k^2 u \end{equation*}

with

\begin{equation*} u(0) = u(L) = 0 \end{equation*}

Following are the codes.

Python

from __future__ import division, print_function
import numpy as np
from scipy.linalg import eigh as eigsh
import matplotlib.pyplot as plt


def beam_FDM_eigs(L, N):
    x = np.linspace(0, L, N)
    dx = x[1] - x[0]
    stiff = np.diag(6*np.ones(N - 2)) -\
            np.diag(4*np.ones(N - 3), -1) - np.diag(4*np.ones(N - 3), 1) +\
            np.diag(1*np.ones(N - 4), 2) + np.diag(1*np.ones(N - 4), -2)
    vals, vecs = eigsh(stiff/dx**4)
    return vals, vecs, x


N = 1001
nvals = 20
nvecs = 4
vals, vecs, x = beam_FDM_eigs(1.0, N)

#%% Plotting
num = np.linspace(1, nvals, nvals)
plt.rcParams["mathtext.fontset"] = "cm"
plt.figure(figsize=(8, 3))
plt.subplot(1, 2, 1)
plt.plot(num, np.sqrt(vals[0:nvals]), "o")
plt.xlabel(r"$N$")
plt.ylabel(r"$\omega\sqrt{\frac{\lambda}{EI}}$")
plt.subplot(1, 2 ,2)
for k in range(nvecs):
    vec = np.zeros(N)
    vec[1:-1] = vecs[:, k]
    plt.plot(x, vec, label=r'$n=%i$'%(k+1))

plt.xlabel(r"$x$")
plt.ylabel(r"$w$")
plt.legend(ncol=2, framealpha=0.8, loc=1)
plt.tight_layout()
plt.show()

Julia

using PyPlot


function beam_FDM_eigs(L, N)
    x = linspace(0, L, N)
    dx = x[2] - x[1]
    stiff = diagm(6*ones(N - 2)) -
            diagm(4*ones(N - 3), -1) - diagm(4*ones(N - 3), 1) +
            diagm(1*ones(N - 4), 2) + diagm(1*ones(N - 4), -2)
    vals, vecs = eig(stiff/dx^4)
    return vals, vecs, x
end


N = 1001
nvals = 20
nvecs = 4
vals, vecs, x = beam_FDM_eigs(1.0, N)

#%% Plotting
num = 1:nvals
# Missing line for setting the math font
figure(figsize=(8, 3))
subplot(1, 2, 1)
plot(num, sqrt.(vals[1:nvals]), "o")
xlabel(L"$N$")
ylabel(L"$\omega\sqrt{\frac{\lambda}{EI}}$")
subplot(1, 2 ,2)
for k in 1:nvecs
    vec = zeros(N)
    vec[2:end-1] = vecs[:, k]
    plot(x, vec, label="n=$(k)")
end

xlabel(L"$x$")
ylabel(L"$w$")
legend(ncol=2, framealpha=0.8, loc=1)
tight_layout()
show()

Both have (almost) the same result, as follows

Vibration modes of an encastred beam.

Comparison Python/Julia

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

For Python:

%timeit beam_FDM_eigs(1.0, N)

with result

10 loops, best of 3: 196 ms per loop

For Julia:

@benchmark beam_FDM_eigs(1.0, N)

with result

BenchmarkTools.Trial:
  memory estimate:  99.42 MiB
  allocs estimate:  55
  --------------
  minimum time:     665.152 ms (17.14% GC)
  median time:      775.441 ms (21.76% GC)
  mean time:        853.401 ms (16.86% GC)
  maximum time:     1.236 s (15.68% GC)
  --------------
  samples:          6
  evals/sample:     1

In this case, we can say that the Python code is roughly 4 times faster than Julia.

Comments

Comments powered by Disqus