In [2]:
from IPython.html.services.config import ConfigManager
from IPython.utils.path import locate_profile
cm = ConfigManager(profile_dir=locate_profile(get_ipython().profile))
cm.update('livereveal', {
'theme': 'sky',
'transition': 'zoom',
'start_slideshow_at': 'selected',
})
Out[2]:
In [9]:
%matplotlib inline
import matplotlib.pyplot as plt
plt.xkcd()
import seaborn as sns
The time-dependent wave equation
$$\frac{1}{c^2} \frac{\partial^2 p}{\partial t^2} = \Delta p + \mbox{boundary conditions}.$$$c$ is the wave speed.
A typical way is to solve for the individual frequency:
$$p(x, t) = \int e^{i w t} p(x) dw,$$giving the scalar Helmholtz problem
$$ \Delta p + k^2 p = 0, +\mbox{boundary conditions}, \quad k = \frac{w}{c}.$$Physical meaning: harmonic source (1 Ghz), and it is often sufficient.
Wideband problems exist, but are more difficult.
It is not easy to create "free space" conditions experimentally. There is a special Anechoic chamber that models such free space.
Using Dirichlet/Neumann we can get resonances (i.e., non-invertible matrices)
In [47]:
import numpy as np
import scipy.linalg
import tt
d = 8
n = 2 ** d
t = np.linspace(0, 1, n + 1) #Left point ommited
h = t[1] - t[0]
a = np.zeros((n, n), dtype=np.complex)
K = np.sqrt(9.79282562485+0j) #For the resonance frequency
#K = 22
M = np.zeros((n, n))
for i in xrange(n-1):
a[i, i] = 2.0/h**2 - K ** 2
a[i, i+1] = -1.0 * (1/h**2)
a[i+1, i] = -1.0 * (1/h**2)
a[n-1, n-1] = 2.0/h**2 - K ** 2#Dirichlet boundary conditions
#rhs = np.random.randn(n)
rhs = np.zeros(n)
#rhs[0] = 1.0
rhs = np.ones(n)
sol = np.linalg.solve(a, rhs)
#print np.min(np.linalg.eigvals(a))
plt.plot(sol.real)
Out[47]:
The right boundary condition in this case is the so-called Robin boundary condition of the form
$$u'(1) + i k \beta u(1) = 0$$
In [49]:
import numpy as np
import scipy.linalg
import tt
d = 8
n = 2 ** d
t = np.linspace(0, 1, n + 1) #Left point ommited
h = t[1] - t[0]
a = np.zeros((n, n), dtype=np.complex)
#K = np.sqrt(9.79282562485+0j) #For the resonance frequency
K = 22
M = np.zeros((n, n))
for i in xrange(n-1):
a[i, i] = 2.0/h**2 - K ** 2
a[i, i+1] = -1.0 * (1/h**2)
a[i+1, i] = -1.0 * (1/h**2)
a[n-1, n-1] = (1/h**2 - K ** 2) * 1.0 - 1j * K
rhs = np.zeros(n)
rhs[0] = 1.0
sol = np.linalg.solve(a, rhs)
#print np.min(np.linalg.eigvals(a))
plt.plot(sol.real)
Out[49]:
Fundamental solution for the Helmholtz equation has the form
$$\Phi_k(x) = \begin{cases} \frac{i}{4} H^{(1)}_0 (k \Vert x - y \Vert), \quad d = 2, \\ \frac{e^{ik \Vert x - y \Vert}}{4 \pi \Vert x - y \Vert}, \quad d = 3, \end{cases} $$and $H^{(1)}_{\nu}$ is the Hankel function of the first kind of order $\mu$
In the case of sound-soft scattering (i.e., Dirichlet boundary conditions of the form)
$$p_{total} = p_{ext} + p_{scat} = 0.$$It is a good idea to use double layer potential (i.e. the unknowns are $\frac{\partial u}{\partial n})$.
$$\frac{1}{2} \frac{\partial u}{\partial n} + \int_{\Gamma} \Big( \frac{\partial \Phi_k(x, y)}{\partial n(x)} - i \eta \Phi_k(x, y) \Big)\frac{\partial u}{\partial n}(y)ds(y) = f_{k, \eta},$$where
$$ f_{k, n} = \frac{\partial u^I}{\partial n}(x) - i \eta u^{I}(x). $$In the operator form we have
$$\Big(\frac{1}{2} I + D'_k - i \eta S_k\Big) v = f_{k, \eta}.$$This is called combined potential equation, and $\eta$ is the parameter.
Note, that due to the complex part, the matrix will be made always non-singular (which is not
the case in the "resonant" frequencies).
If solved by standard methods, the matrices will be dense, and the number of unknowns will grow with respect to $k$.
Also, the ranks will grow with respect to $k$, at least as $\mathcal{O}(\sqrt{k})$ in the 3D case, so the H-matrix argument will fail.
We need something special to get rid of $k$-dependence!
For the case of very large $k$ we get the case of geometric optics: we either have "shadow region" or do not have it.
However, the most interesting happens, when there is something in the "shadow region", i.e. diffraction.
However, we may use geometric optics to get the right expansion.
One of the possible ways is to use Kirchoff approximation (physical optics) Consider a plane wave $u^I(x) = -\exp(ikx \cdot \widehat{a}$ incident on an infinite plane with unit norm $n$. assumption of sound-soft scattering, the reflected plane wave is
$$ u^R(x) = -\exp(ikx \widehat{a}^R), \quad \widehat{a}^R = \widehat{a} - 2(n, a) n. $$The normal derivative is then
$$ v(x) = 2 \frac{\partial u^I}{\partial n} = 2ikn \widehat{a} \exp{ikx \widehat{a}}. $$Then we use it in the illuminated region, and set $v$ to zero in the shadow region.
Besides the scalar Helmholtz equations, there are the Maxwell equations that describe electromagnetics.
$\nabla \cdot D = 4 \pi \rho$
$\nabla \cdot B = 0$
$\nabla \times E = -\frac{1}{c}\frac{\partial B}{\partial t}$
$\nabla \times H = \frac{1}{c}\Big(4 \pi J + \frac{\partial D}{\partial t}\Big)$
$D$ and $H$ are auxiliary fields in the macroscopic Maxwell equations $D = \varepsilon E, \quad H = \frac{1}{\mu} B$ constitutive relations
We get $\nabla \times \Big(E + \frac{1}{c} i w \mu A\Big) = 0$,
$E + \frac{1}{c} i w \mu A = -\nabla \Phi$
$\nabla \times \nabla \times A - k^2 A = \frac{1}{c} 4 \pi J - \frac{1}{c} i w \varepsilon \nabla \Phi$
and the Lorentz condition
$\nabla \cdot A = - \frac{1}{c} i w \varepsilon \nabla \Phi$
$\Delta A + k^2 A = -\frac{4 \pi}{c} J$
That is the starting point for the Electric Field Integral equation.
In [4]:
from IPython.core.display import HTML
def css_styling():
styles = open("./styles/custom.css", "r").read()
return HTML(styles)
css_styling()
Out[4]: