In [337]:
using PyPlot
PyPlot.version
Out[337]:
For now, we want to solve the diffusion problem on a rectangular domain. The reaction part of the equation we want to solve in the end is currently set to zero. We therefore want to find a function $u$ defined as follows.
\begin{align} u: [a, b] \times [c, d] \times ℝ^+ &\to ℝ^+ \\ (x, y, t) &\to u(x, y, t) \end{align}We therefore want to solve the following PDE:
$$ u_t = D Δ u = D (u_xx + u_yy) \\ u(a,y,t) = u(b,y,t) = u(x,c,t)=u(x,d,t)=0 $$For this first solution, we use forward Euler integration in time and a second-order finite differences method in space. This leads to the following method to calculate the next time step:
$$ u^{n+1} = u^n + D \cdot \Delta t \cdot A u^n $$Here, $A$ is the finite differences matrix for the 2D Laplace operator.
In [338]:
D = 0.01 # diffusion coefficient;
In [339]:
# define domain of the problem (rectangle) and time interval
xstart = 1
xend = 3
ystart = 2
yend = 5
tstart = 0
tend = 10
Out[339]:
In [340]:
# define number of intervals in space
Nx = 64
Ny = 64
Out[340]:
In [341]:
dx = (xend - xstart) / Nx
dy = (yend - ystart) / Ny
Out[341]:
In order to have a stable solution, we set $\Delta t$ based on $\Delta x$ and $\Delta y$, as follows:
$$ \Delta t \propto \frac{\Delta x^2 + \Delta y^2}{D} $$where $D$ is the diffusion coefficient.
In [342]:
securityfactor = 10 # at least 9 is needed (with Nx=Ny=64)
Nt = int(securityfactor * (tend - tstart) * D / (dx^2 + dy^2))
dt = (tend - tstart) / Nt
Out[342]:
As a very simple scenario, we use the following initial conditions:
$$ u(x,y,t = 0) = C_0 \cdot sin(π \frac{x-a}{b-a}) \cdot sin(π \frac{y-c}{d-c}) $$This leads to the following analytical solution:
$$ u(x,y,t) = C_0 \cdot sin( \frac{x-a}{b-a}) \cdot sin(π \frac{y-c}{d-c}) \cdot exp( -\left( \frac{1}{(b-a)^2} + \frac{1}{(d-c)^2} \right) \pi^2 D t) $$
In [343]:
C0 = 5
exact(x,y,t) = C0 * sin((x-xstart)/(xend-xstart) * pi) * sin((y-ystart)/(yend-ystart) * pi) *
exp( - (1/(xend-xstart)^2 + 1/(yend-ystart)^2) * pi^2 * D * t)
Out[343]:
In [344]:
C_exact = Float64[exact(x,y,tend) for x=xstart:dx:xend, y=ystart:dy:yend]
pcolormesh(C_exact)
colorbar()
clim((0,C0))
In [345]:
# use exact function to set up initial conditions
C = Float64[exact(x,y,0) for x=xstart:dx:xend, y=ystart:dy:yend];
We define the finite differences matrices seperately for the $\partial_{xx}$ and $\partial_{yy}$ operators. The vector of unknowns $u$ is organized in col-major fashion (neighboring elements have the same y-value but not the same x-value).
Thus, the matrix for $\partial_{xx}$ is block-diagonal and organized as follows:
\begin{align*} FD_x = \begin{pmatrix} B_1 & 0 & \cdots & 0 \\ 0 & B_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & B_{N_y} \end{pmatrix} & \hspace{1cm} \text{with} & B_i = \begin{pmatrix} -2 & 1 & 0 & \cdots & 0 \\ 1 & -2 & 1 & \ddots & \vdots \\ 0 & 1 & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & 1 \\ 0 & \cdots & 0 & 1 & -2 \end{pmatrix}_{N_x \times N_x} \end{align*}The matrix for $\partial_{xx}$ has has non-zero entries in the diagonal and the $N_{xx}$-th upper and lower diagonals.
\begin{align*} FD_y = \begin{pmatrix} -2 & 0 & \cdots & 1 & \cdots & 0 \\ 0 & -2 & 0 & \ddots & \ddots & \vdots \\ \vdots & 0 & -2 & \ddots & \ddots & 1 \\ 1 & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & 1 & \cdots & 0 & -2 \end{pmatrix}_{N_x N_y \times N_x N_y} \end{align*}The sum of these matrices, combined with their step sizes, forms the finite difference matrix that is used for the integration.
$$ FD = \frac{1}{\Delta x^2} FD_x + \frac{1}{\Delta y^2} FD_y $$
In [346]:
# number of unknowns in x and y direction
ux = (Nx-1);
uy = (Ny-1);
In [347]:
# finite difference matrix for ∂u/∂x
FDx_local = spdiagm((ones(ux-1),-2*ones(ux),ones(ux-1)),(-1,0,1));
FDx = blkdiag({FDx_local for i in 1:uy}...);
# finite difference matrix for ∂u/∂y
FDy = spdiagm((ones(ux*(uy-1)),-2*ones(ux*uy),ones(ux*(uy-1))),(-ux,0,ux));
In [348]:
FD = 1/dx^2 * FDx + 1/dy^2 * FDy;
In [349]:
u = vcat(C[2:end-1,2:end-1]...);
iterationmatrix = dt * D * FD;
In [350]:
for i in 1:Nt
u += iterationmatrix * u
end
C[2:end-1,2:end-1] = reshape(u,ux,uy)
pcolormesh(C)
colorbar()
clim((0,C0))
In [351]:
vecnorm(C - C_exact)
Out[351]:
In [352]:
# collected values for varying x (Ny = 64, Nt = 10000)
error_x = [
0.3506375511215073 # Nx = 4
0.12577467828850733 # Nx = 8
0.045244709643395255 # Nx = 16
0.016937758467202517 # Nx = 32
0.0073057790525524555 # Nx = 64
0.004444809104674145 # Nx = 128
]
Nxs = [4 8 16 32 64 128]'
loglog(Nxs, error_x)
grid(true)