```
```{raw} latex
%%start demo%%
```
```{code-cell}
A = [1 2; -2 0]
```
```{code-cell}
B = [ 1 10 100; -5 5 3 ]
```
Applying the definition manually, we get
```{code-cell}
A_kron_B = [ A[1,1]*B A[1,2]*B;
A[2,1]*B A[2,2]*B ]
```
```{index} ! Julia; kron
```
But `kron` is a built-in function.
```{code-cell}
kron(A,B)
```
```{raw} html

```
```{raw} latex
%%end demo%%
```
The Kronecker product obeys several natural-looking identities:
(theorem-laplace-kron)=
::::{proof:theorem} Kronecker product identities
Given matrices for which the non-Kronecker operations make sense, the following hold.
1. $\mathbf{A}\otimes (\mathbf{B} + \mathbf{C}) = \mathbf{A}\otimes \mathbf{B} + \mathbf{A}\otimes \mathbf{C}$
2. $(\mathbf{A} + \mathbf{B}) \otimes \mathbf{C} = \mathbf{A}\otimes \mathbf{C} + \mathbf{B}\otimes \mathbf{C}$
3. $(\mathbf{A} \otimes \mathbf{B}) \otimes \mathbf{C} = \mathbf{A} \otimes (\mathbf{B} \otimes \mathbf{C})$
4. $(\mathbf{A} \otimes \mathbf{B})^T = \mathbf{A}^T \otimes \mathbf{B}^T$
5. $(\mathbf{A} \otimes \mathbf{B})^{-1} = \mathbf{A}^{-1} \otimes \mathbf{B}^{-1}$
6. $(\mathbf{A} \otimes \mathbf{B})(\mathbf{C}\otimes \mathbf{D}) = (\mathbf{A}\mathbf{C}) \otimes (\mathbf{B}\mathbf{D})$
For the vec operation defined in {numref}`Definition {number}
```
```{raw} latex
%%start demo%%
```
Here is a forcing function for Poisson's equation.
```{code-cell}
f = (x,y) -> x^2 - y + 2;
```
We make a crude discretization for illustrative purposes.
```{code-cell}
m,n = 6,5
x,Dx,Dxx = FNC.diffmat2(m,[0,3])
y,Dy,Dyy = FNC.diffmat2(n,[-1,1])
unvec = u -> reshape(u,m+1,n+1);
```
Next, we evaluate $\phi$ on the grid.
```{code-cell}
F = [ f(x,y) for x in x, y in y ]
```
Here are the equations for the PDE collocation, before any modifications are made for the boundary conditions.
```{code-cell}
A = kron(I(n+1),Dxx) + kron(Dyy,I(m+1))
b = vec(F);
```
The number of equations is equal to $(m+1)(n+1)$, which is the total number of points on the grid.
```{code-cell}
@show N = length(F);
```
The combination of Kronecker products and finite differences produces a characteristic sparsity pattern.
```{code-cell}
spy(sparse(A),color=:blues,m=3,
title="System matrix before boundary conditions")
```
We now construct a Boolean array the same size as `F` to indicate where the boundary points lie in the grid.
```{code-cell}
isboundary = trues(m+1,n+1)
isboundary[2:m,2:n] .= false
idx = vec(isboundary);
```
```{code-cell}
spy(sparse(isboundary),m=3,color=:darkblue,legend=:none,
title="Boundary points",
xaxis=("column index",[0,n+2]),yaxis=("row index",[0,m+2]) )
```
In order to impose Dirichlet boundary conditions, we replace the boundary rows of the system by rows of the identity.
```{code-cell}
I_N = I(N)
A[idx,:] .= I_N[idx,:]; # Dirichlet conditions
```
```{code-cell}
spy(sparse(A),color=:blues,m=3,
title="System matrix with boundary conditions")
```
Finally, we must replace the rows in the vector $\mathbf{b}$ by the boundary values being assigned to the boundary points. Here, we let the boundary values be zero everywhere.
```{code-cell}
b[idx] .= 0; # Dirichlet values
```
Now we can solve for $\mathbf{u}$ and reinterpret it as the matrix-shaped $\mathbf{U}$, the solution on our grid.
```{code-cell}
u = A\b
U = unvec(u)
```
```{raw} html

```
```{raw} latex
%%end demo%%
```
## Implementation
{numref}`Function {number}
```
```{raw} latex
%%start demo%%
```
We can engineer an example by choosing the solution first. Let $u(x,y)=\sin(3xy-4y)$. Then one can derive $f=\Delta u = -\sin(3xy-4y)\bigl(9y^2+(3x-4)^2\bigr)$ for the forcing function and use $g=u$ on the boundary.
First we define the problem on $[0,1]\times[0,2]$.
```{code-cell}
f = (x,y) -> -sin(3*x*y-4*y)*(9*y^2+(3*x-4)^2);
g = (x,y) -> sin(3*x*y-4*y);
xspan = [0,1]; yspan = [0,2];
```
Here is the finite-difference solution.
```{code-cell}
x,y,U = FNC.poissonfd(f,g,40,xspan,60,yspan);
```
```{code-cell}
:tags: [hide-input]
surface(x,y,U',color=:viridis,
title="Solution of Poisson's equation",
xaxis=(L"x"),yaxis=(L"y"),zaxis=(L"u(x,y)"),
right_margin=3Plots.mm,camera=(70,50))
```
The error is a smooth function of $x$ and $y$. It must be zero on the boundary; otherwise, we have implemented boundary conditions incorrectly.
```{code-cell}
:tags: [hide-input]
error = [g(x,y) for x in x, y in y] - U;
M = maximum(abs,error)
contour(x,y,error',levels=17,aspect_ratio=1,
clims=(-M,M),color=:redsblues,colorbar=:bottom,
title="Error",xaxis=(L"x"),yaxis=(L"y"),
right_margin=7Plots.mm)
plot!([0,1,1,0,0],[0,0,2,2,0],l=(2,:black))
```
```{raw} html

```
```{raw} latex
%%end demo%%
```
## Accuracy and efficiency
In {numref}`Function {number}