Ir al contenido principal

# Numerical methods challenge: Day 26

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.

## Boundary element method

Today we have the Boundary element method , or, an attempt. We want so solve the equation

\begin{equation*} \nabla^2 u = -f(x) \end{equation*}

with

\begin{equation*} u(x) = g(x),\quad \forall (x, y)\in \partial \Omega \end{equation*}

For this method, we need to use an integral representation of the equation, that in this case is

\begin{equation*} u(\xi) = \int_{S} [u(x) F(x, \xi) - q(x)G(x, \xi)]dS(x) + \int_{V} f(x) G(x, \xi) dV(x) \end{equation*}

with

\begin{equation*} G(x,\xi)= -\frac{1}{2\pi}\ln|x- \xi| \end{equation*}

and

\begin{equation*} F(x,\xi) = -\frac{1}{2\pi |x- \xi|^2}(x - \xi)\cdot\hat{n} \end{equation*}

Then, we can form a system of equations

\begin{equation*} [G]\{q\} = [F]\{u\}\, , \end{equation*}

that we obtain by discretization of the boundary. If we take constant variables over the discretization, the integral can be computed analytically by

\begin{equation*} G_{nm} = -\frac{1}{2\pi}\left[r \sin\theta\left(\ln|r| - 1\right) + \theta r\cos\theta\right]^{\theta_B, r_B}_{\theta_A, r_A} \end{equation*}

and

\begin{equation*} F_{nm} = \left[\frac{\theta}{2\pi}\right]^{\theta_B}_{\theta_A} \end{equation*}

for points $n$ and $m$ in different elements, and the subindices $A,B$ refer to the endpoints of the evaluation element. For diagonal terms the integrals evaluate to

\begin{equation*} G_{nn} = -\frac{L}{2\pi}\left(\ln\left\vert\frac{L}{2}\right\vert - 1\right) \end{equation*}

and

\begin{equation*} F_{nn} = - \frac{1}{2\pi} \end{equation*}

with $L$ the size of the element.

Below is the Python code. In this case, I did not succeed and the code is not working right now.

from __future__ import division, print_function
import numpy as np
from numpy import log, sin, cos, arctan2, pi, mean, arctan
from numpy.linalg import norm, solve
import matplotlib.pyplot as plt

def assem(coords, elems):
nelems = elems.shape
Gmat = np.zeros((nelems, nelems))
Fmat = np.zeros((nelems, nelems))
for ev_cont, elem1 in enumerate(elems):
print(coords[elem1], coords[elem1])
for col_cont, elem2 in enumerate(elems):
pt_col = mean(coords[elem2], axis=0)
if ev_cont == col_cont:
L = norm(coords[elem1] - coords[elem1])
Gmat[ev_cont, ev_cont] = - L/(2*pi)*(log(L/2) - 1)
Fmat[ev_cont, ev_cont] = - 0.5
else:
Gij, Fij = influence_coeff(elem1, coords, pt_col)
Gmat[ev_cont, col_cont] = Gij
Fmat[ev_cont, col_cont] = Fij
return Gmat, Fmat

def influence_coeff(elem, coords, pt_col):
r_A = coords[elem] - pt_col
r_B = coords[elem] - pt_col
theta_A = arctan2(r_A, r_A)
theta_B = arctan2(r_B, r_B)
G_coeff = r_B*(log(norm(r_B)) - 1) + theta_B*r_B -\
(r_A*(log(norm(r_A)) - 1) + theta_A*r_A)
F_coeff = theta_B - theta_A
return -G_coeff/(2*pi), F_coeff/(2*pi)

def eval_sol(ev_coords, coords):
nelems = elems.shape
Gmat = np.zeros((nelems, nelems))
Fmat = np.zeros((nelems, nelems))
for ev_cont, elem1 in enumerate(elems):
L = norm(coords[elem1] - coords[elem1])
for col_cont, elem2 in enumerate(elems):
pt_col = mean(coords[elem2], axis=0)
if ev_cont == col_cont:
Gmat[ev_cont, ev_cont] = - L/(2*pi)*(log(L/2) - 1)
Fmat[ev_cont, ev_cont] = - 0.5
else:
Gmat[ev_cont, col_cont], Fmat[ev_cont, col_cont] = \
influence_coeff(elem1, coords, pt_col)

nelems = 3