Ir al contenido principal

# Numerical methods challenge: Day 22

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

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

## Comentarios

Comments powered by Disqus
{"extraction_info": {"found_math": true, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 1, "mathjax_display_tex": 0, "mathjax_asciimath": 1, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0, "katex": 0, "math-container": 0, "wp-katex-eq": 0, "align": 0, "equation": 0, "x-ck12": 0, "texerror": 0, "math_score": 0.2611284852027893, "perplexity": 20687.44781548873}, "config": {"markdown_headings": true, "markdown_code": true, "boilerplate_config": {"ratio_threshold": 0.18, "absolute_threshold": 10, "end_threshold": 15, "enable": false}, "remove_buttons": true, "remove_image_figures": true, "remove_link_clusters": true, "table_config": {"min_rows": 2, "min_cols": 3, "format": "plain"}, "remove_chinese": true, "remove_edit_buttons": true, "extract_latex": true}, "warc_path": "s3://commoncrawl/crawl-data/CC-MAIN-2020-40/segments/1600400189928.2/warc/CC-MAIN-20200919013135-20200919043135-00048.warc.gz"}