Skip to main content

Nicolás' blog

Numerical methods challenge: Day 21

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.

Jacobi iteration

Today we have a Finite difference method combined with Jacobi method We are solving the Laplace equation.

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

with

\begin{align*} u(0, y) = 1 -y,\quad u(x, 0) = 1 - x,\\ u(1, y) = u(x, 1) = 0 \end{align*}

Following are the codes.

Python

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


def jacobi(u, update, tol=1e-6, niter=500):
    for n in range(niter):
        u_new = update(u)
        if np.linalg.norm(u_new - u) < tol:
            break
        else:
            u = u_new.copy()
    return u, n


def heat_FDM(u):
    u_new = u.copy()
    u_new[1:-1, 1:-1] = 0.25*(u_new[0:-2, 1:-1] + u_new[2:, 1:-1] +\
         u_new[1:-1, 0:-2] + u_new[1:-1, 2:])
    return u_new


nx = 50
ny = 50
x_vec = np.linspace(-0.5, 0.5, nx)
y_vec = np.linspace(-0.5, 0.5, ny)
u0 = np.zeros((nx, ny))
u0[:, 0] = 1 - x_vec
u0[0, :] = 1 - y_vec
nvec = [100, 1000, 10000, 100000]
for num, niter in enumerate(nvec):
    u, n = jacobi(u0, heat_FDM, tol=1e-12, niter=niter)
    plt.subplot(2, 2, num + 1)
    plt.contourf(u, cmap='hot')
    plt.xlabel(r"$x$")
    plt.ylabel(r"$y$")
    plt.title('{} iterations'.format(n + 1))
    plt.axis('image')
plt.tight_layout()
plt.show()

Julia

using PyPlot


function jacobi(u, update; tol=1e-6, niter=500)
    num = niter
    for n = 1:niter
        u_new = update(u)
        if norm(u_new - u) < tol
            num = n
            break
        else
            u = copy(u_new)
        end
    end
    return u, num
end


function heat_FDM(u)
    u_new = copy(u)
    u_new[2:end-1, 2:end-1] = 0.25*(u_new[1:end-2, 2:end-1] +
        u_new[3:end, 2:end-1] + u_new[2:end-1, 1:end-2] + u_new[2:end-1, 3:end])
    return u_new
end


nx = 50
ny = 50
x_vec = linspace(-0.5, 0.5, nx)
y_vec = linspace(-0.5, 0.5, ny)
u0 = zeros(nx, ny)
u0[:, 1] = 1 - x_vec
u0[1, :] = 1 - y_vec
nvec = [100, 1000, 10000, 100000]
for (num, niter) = enumerate(nvec)
    u, n = jacobi(u0, heat_FDM, tol=1e-12, niter=niter)
    subplot(2, 2, num)
    contourf(u, cmap="hot")
    xlabel(L"$x$")
    ylabel(L"$y$")
    title("$(n) iterations")
    axis("image")
end
tight_layout()
show()
Solution of the differential equation that satisfy the boundary conditions.

Comparison Python/Julia

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

For Python:

%timeit jacobi(u0, heat_FDM, tol=1e-12, niter=1000)

with result

10 loops, best of 3: 29.6 ms per loop

For Julia:

@benchmark jacobi(u0, heat_FDM, tol=1e-12, niter=1000)

with result

BenchmarkTools.Trial:
  memory estimate:  247.89 MiB
  allocs estimate:  43002
  --------------
  minimum time:     196.943 ms (5.66% GC)
  median time:      203.230 ms (5.74% GC)
  mean time:        203.060 ms (6.01% GC)
  maximum time:     208.017 ms (5.51% GC)
  --------------
  samples:          25
  evals/sample:     1

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

Comments

Comments powered by Disqus