Elliptic Partial Differential Equations

3. Elliptic Partial Differential Equations#

import matplotlib.pyplot as plt
import numpy as np
from helper import Approx
from scipy.sparse import spdiags

# from scipy.sparse.linalg import spsolve
from sympy import Derivative as D
from sympy import Eq, Function, Symbol, solve, symbols

x, y = symbols("x, y")
U = symbols("U", cls=Function)(x, y)
f = symbols("f", cls=Function)(x, y)

3.1. Laplace’s equation#

Eq(D(U, x, 2) + D(U, y, 2), 0)
\[\displaystyle \frac{\partial^{2}}{\partial x^{2}} U{\left(x,y \right)} + \frac{\partial^{2}}{\partial y^{2}} U{\left(x,y \right)} = 0\]

3.2. Poisson’s equation#

Eq(D(U, x, 2) + D(U, y, 2), f)
\[\displaystyle \frac{\partial^{2}}{\partial x^{2}} U{\left(x,y \right)} + \frac{\partial^{2}}{\partial y^{2}} U{\left(x,y \right)} = f{\left(x,y \right)}\]

3.3. Discretising a two-dimensional domain#

\[ \Delta x= \dfrac{x_{\max}-x_{\min}}{N_x-1} \]
\[ \Delta y= \dfrac{y_{\max}-y_{\min}}{N_y-1} \]
\[ x_{i}= x_{\min}+i\Delta x \]
\[ y_{j}= y_{\min}+j\Delta y \]
fig = plt.figure()
for i in range(3):
    plt.plot((0, 2), (i, i), "k")

for j in range(3):
    plt.plot((j, j), (0, 2), "k")

plt.plot(0, 0, "ko")
plt.plot(1, 0, "ko")
plt.plot(2, 0, "ko")
plt.plot(0, 1, "ko")
plt.plot(1, 1, "ko")
plt.plot(2, 1, "ko")
plt.plot(0, 2, "ko")
plt.plot(1, 2, "ko")
plt.plot(2, 2, "ko")
ax = fig.gca()
ax.axis("off")
ax.axis("equal")
plt.text(0, -0.15, "$i-1$", ha="center", fontsize=12)
plt.text(1, -0.15, "$i$", ha="center", fontsize=12)
plt.text(2, -0.15, "$i+1$", ha="center", fontsize=12)
plt.text(-0.25, 0, "$j-1$", ha="center", fontsize=12)
plt.text(-0.25, 1, "$j$", ha="center", fontsize=12)
plt.text(-0.25, 2, "$j+1$", ha="center", fontsize=12)
plt.text(0.05, 0.1, "$u^{j-1}_{i-1}$", fontsize=14)
plt.text(1.05, 0.1, "$u^{j-1}_{i}$", fontsize=14)
plt.text(2.05, 0.1, "$u^{j-1}_{i+1}$", fontsize=14)
plt.text(0.05, 1.1, "$u^{j}_{i-1}$", fontsize=14)
plt.text(1.05, 1.1, "$u^{j}_{i}$", fontsize=14)
plt.text(2.05, 1.1, "$u^{j}_{i+1}$", fontsize=14)
plt.text(0.05, 2.1, "$u^{j+1}_{i-1}$", fontsize=14)
plt.text(1.05, 2.1, "$u^{j+1}_{i}$", fontsize=14)
plt.text(2.05, 2.1, "$u^{j+1}_{i+1}$", fontsize=14);
# plt.text(-0.45, 2, "$y$")
# plt.text(2.5, 0, "$x$")
../_images/c753cfb12a773a6f107308992182462032fe216f704efb26d1703e2850e78d05.png

3.3.1. Boundary conditions#

3.3.1.1. Dirichlet boundary conditions#

\[\begin{split} \begin{aligned} u_{0}&=\alpha\\ u_{N-1}&=\beta \end{aligned} \end{split}\]
\[\begin{split} \begin{aligned} u_{0}&=f\left(t,x,y\right)\\ u_{N-1}&=g\left(t,x,y\right) \end{aligned} \end{split}\]
\[\begin{split} \begin{aligned} u_{0}&=u_{N-2}\\ u_{N-1}&=u_{1} \end{aligned} \end{split}\]

3.3.1.2. Neumann boundary conditions#

\[\begin{split} \begin{aligned} u_{-1}&=u_{1}\\ u_{N}&=u_{N-2} \end{aligned} \end{split}\]
\[\begin{split} \begin{aligned} v_{-1}&=-v_{1}\\ v_{N}&=-v_{N-2} \end{aligned} \end{split}\]
f = symbols("f", cls=Function)(U)

Eq(D(U, x), f)
\[\displaystyle \frac{\partial}{\partial x} U{\left(x,y \right)} = f{\left(U{\left(x,y \right)} \right)}\]

3.3.2. Deriving a FDS to solve Laplace’s equation#

x, y = symbols("x y")
dx = Symbol(r"\Delta x")
dy = Symbol(r"\Delta y")
U = Function("U")

_ = U(x + dx, y).series(x=dx, x0=0, n=4).removeO().simplify()
__ = U(x - dx, y).series(x=dx, x0=0, n=4).removeO().simplify()
U_xx = Approx(
    D(U(x, y), x, 2),
    solve(Eq(_ + __, U(x + dx, y) + U(x - dx, y)), D(U(x, y), x, 2))[0],
)

U_xx
\[\displaystyle \frac{\partial^{2}}{\partial x^{2}} U{\left(x,y \right)} \approx \frac{- 2 U{\left(x,y \right)} + U{\left(- \Delta x + x,y \right)} + U{\left(\Delta x + x,y \right)}}{\Delta x^{2}}\]
_ = U(x, y + dy).series(x=dy, x0=0, n=4).removeO().simplify()
__ = U(x, y - dy).series(x=dy, x0=0, n=4).removeO().simplify()
U_yy = Approx(
    D(U(x, y), y, 2),
    solve(Eq(_ + __, U(x, y + dy) + U(x, y - dy)), D(U(x, y), y, 2))[0],
)

U_yy
\[\displaystyle \frac{\partial^{2}}{\partial y^{2}} U{\left(x,y \right)} \approx \frac{- 2 U{\left(x,y \right)} + U{\left(x,- \Delta y + y \right)} + U{\left(x,\Delta y + y \right)}}{\Delta y^{2}}\]

3.3.2.1. Second-order finite-difference scheme to solve Laplace’s equation#

\[ u^{j}_{i}= \dfrac{{\left(\Delta y\right)}^{2}\left(u^{j}_{i-1}+u^{j}_{i+1}\right)+{\left(\Delta x\right)}^{2}\left(u^{j-1}_{i}+u^{j+1}_{i}\right)} {2\left({\left(\Delta x\right)}^{2}+{\left(\Delta y\right)}^{2}\right)} \]

3.3.2.2. Solving Laplace’s equation#

\[ 0\leq x\leq 4,\quad 0\leq y\leq 3. \]
N_x, N_y = 5, 4
x, dx = np.linspace(start=0, stop=4, num=N_x, retstep=True)
y, dy = np.linspace(start=0, stop=3, num=N_y, retstep=True)

x

array([0., 1., 2., 3., 4.])

y

array([0., 1., 2., 3.])

def Creating_Matrix_A(Nx, Ny, r):
    # Diagonal principal
    diag0 = np.tile(np.full(Nx, 1 - 4 * r), Nx)
    # Diagonal principal superior
    diag1 = np.tile(np.concatenate((np.array([0.0, 2 * r]), np.full(Nx - 2, r))), Nx)
    # Diagonal principal inferior
    diag_1 = np.flip(diag1)

    # Diagonal superior
    diag_sup = np.concatenate(
        (np.zeros(Nx), np.full(Nx, 2 * r), np.tile(np.full(Nx, r), Nx - 2))
    )
    # Diagonal inferior
    diag_inf = np.flip(diag_sup)

    # Retornando Matriz sparse
    return spdiags(
        [diag_inf, diag_1, diag0, diag1, diag_sup],
        [-Nx, -1, 0, 1, Nx],
        Nx**2,
        Ny**2,
    ).tocsr()

A = Creating_Matrix_A(2, 2, 0.5)

A.todense()

matrix([[-1.,  1.,  1.,  0.],
        [ 1., -1.,  0.,  1.],
        [ 1.,  0., -1.,  1.],
        [ 0.,  1.,  1., -1.]])
Contents