Motivation

Why computer modelling?

Because it is cheaper than real-life experiment, or in the case when real-life experiment is not possible.

Typical steps of computer modelling

  1. Formulate the mathematical problem as an equation for some quantities

  2. Replace the continious problem by a discrete one (discretization)

  3. Solve the resulting discrete problem

The simplest cycle: Mathematical model $\to$ Discretization $\to$ Solve

Discretization

The discretization is replacement of the region by discrete elements:

Random notes

  • Discretization and solve can be connected
  • Fast solvers are needed
  • Only a subproblem in design and optimization
  • Many physical problems are still too complex for math (turbulence!)

Consider

It takes a lot to create

  1. A model
  2. Discretization
  3. Solvers

What if the computer time to compute a forecast for 1 day is more than 1 day?

Many process in physics are modelled as PDEs.

  • Diffusion processes (heat transfer), electrostatics (circuit design) Poisson equation
  • Sound propagation (noise on the streets, buildings) – Helmholtz equation
  • Electromagnetics – MRI (magnetic resonance imaging) – Maxwell equations
  • Fluid flows – Stokes / Navier Stokes equations

These are all partial differential equations!

PDEs appear in many areas, including

  • Modelling of physical processes (heat, elasticity, fluid flows)
  • Financial math (Black Scholes equation)
  • Chemical engineering (Smoluchowsky equation)
  • Nanoworld (Schrodinger equation)
  • Optimal control of robots (Hamilton-Jacobi-Bellman equation)

Why do we need fast methods?

Because the growth of the computer power, being exponential, is still not enough!

The growth due to the algorithm improvements is comparable (and the human brain is still the most energy-efficient computing element)

Source

What do we mean by fast methods?

By fast methods we mean improving the asymptotics with respect to the problem size.

Consider solution of linear system with a sparse matrix $A$:

$$Au = b,$$

where $A$ is a 5-point Laplacian discretization:

$$\frac{u_{i+1, j} + u_{i-1, j} + u_{i, j-1} + u_{i, j+1} - 4 u_{ij} }{h^2} = f_{ij}.$$

What are the complexities (next slide, but let us guess).

Complexity ( essentials )

  • Dense Gaussian elimination: $\mathcal{O}(N^3)$, works up to $10^4$
  • Sparse Gaussian elimination: $\mathcal{O}(N^{\frac{3}{2}})$, works up to $10^6$
  • FFT methods: $\mathcal{O}(N \log N)$, up to $10^8$
  • Multigrid method: $\mathcal{O}(N)$, up to $10^8$
  • Tensor methods for particular right-hand sides (for example, $f=1$), works up to astronomically large sizes ($N = 2^{50}$).

Integral equations

Now, to integral equations!

Physics is described by PDEs

The physics of our world is typically described by local conservation laws, expressed in terms of partial differential equations.

The Poisson equation writes as

$$\Delta u \equiv \text{div} \nabla u = f,$$

plus boundary conditions.

Model problem (electrostatics)

Suppose physical setting: we have an ideally conducting surface $\Omega$ (for example, surface of a cube), which is attached to a battery.

The charges can appear only at the surface, i.e.

$$\Delta V(x) = 0, \quad x\not\in \partial \Omega$$

but at the surface the potential should be constant:

$$V(x) = V_0, \quad x\in \partial\Omega$$
  • This is a classical example of external problem.

  • The potential has to be defined in the full $\mathbb{R}^3$ space.

  • It is quite expensive, boundary conditions on the outer boundary are not straightforward.

From electrostatics to integral equations

The concept of equivalent sources leads to the boundary integral formulation of the problem.

The charges can appear only at the boundary. The charge creates the field $\frac{1}{r}$.

In $\mathbb{R}^3$ function $G(x,y) =\frac{1}{4\pi\|x-y\|}$ is fundamental solution of the operator $\Delta$, since it satisfies

$$\Delta G (x,y) = \delta(x-y),$$

where $\delta$ is a delta-function.

BEM

The boundary integral equation comes from the idea to look for the solution as

$$V(x) = \int_{\partial \Omega} \frac{q(y)}{\Vert x - y\Vert} dy.$$

(it is also called single-layer potential).

Properties

We have $$\Delta V(x) = \int_{\partial \Omega} q(y) \Delta_x\left(\frac{1}{\Vert x - y\Vert}\right) dy = 4\pi\int_{\partial \Omega} q(y) \delta (x-y) dy = (\text{why?}) = 0, \quad x\not\in\partial\Omega$$

therefore it is sufficient to find the unknown function $q$ that satisfies the Dirichlet boundary condition

$$\int_{\partial \Omega} \frac{q(y)}{\Vert x - y\Vert} dy = V_0, \quad x \in \partial \Omega$$

That is the first kind integral equation with singular kernel.

The main benefit is that the unknown function is defined only on $\partial \Omega$!

However, the operator is "non-local" compared with the PDE formulation.

Model problem: acoustics

Room acoustics (for opera houses), noise assessment (for roads, building construction, railways) can be modelled in the same fashion

Model problem: acoustics

The underlying equation is the Helmholtz equation

$$\Delta p + k^2 p = f, $$

plus boundary conditions (typically, Neumann boundary conditions), and $f$ are sound sources (typically, point sources).

The fundamental solution is

$$G(x,y) = \frac{\exp(i k \Vert x - y \Vert)}{\Vert x - y \Vert}.$$

Summary

  • Intro lecture
  • First integral equations

Next lecture

  • How to discretize IE (Nystrom, collocation, Galerkin method, other type of kernels)
  • What are the problems with such discretization.

In [1]:
from IPython.core.display import HTML
def css_styling():
    styles = open("./styles/custom.css", "r").read()
    return HTML(styles)
css_styling()


Out[1]: