`: the boundary conditions are rescaled to read $\sigma u(x,y)=\sigma g(x,y)$, where $\sigma$ is the largest element of a row of $\mathbf{A}$. This tweak improves the condition number of the final matrix.
(demo-laplace-poisson)=
```{proof:demo}
```
```{raw} html
```
```{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} ` we used second-order finite differences for the discretization. Let's simplify the discussion by assuming that $m=n$. The error in the solution can be expected to be $O(n^{-2})$.
The matrix $\mathbf{A}$ has size $N=O(n^2)$. The upper and lower bandwidths of the matrix $\mathbf{A}$ are both $O(n)$. It can be shown that the matrix is symmetric and negative definite, so sparse Cholesky factorization can be applied to $-\mathbf{A}$, taking $O(n^2N)=O(n^4)$ operations.
As $n$ increases, the truncation error decreases while the operation count increases. The growth in operations is faster than the corresponding decrease in error, making it costly to get high accuracy. Suppose we have fixed running time $T$ at our disposal. Then $n=O(T^{1/4})$, and the convergence is $O(1/\sqrt{T})$. For example, reduction of the error by a factor of 10 requires 100 times the computational effort.
If we chose a Chebyshev spectral discretization instead, the calculus changes. Provided that the solution is smooth, we expect convergence at a rate $K^{-n}$ for some $K>1$. However, the system matrix is no longer sparse nor symmetric, and solution by LU factorization takes $O(N^3)=O(n^6)$ flops. Hence as a function of running time $T$, we expect a convergence rate on the order of $K^{-T^{1/6}}$. It's not simple to compare this directly to the finite-difference case. In the long run, the spectral method will be much more efficient, but if the accuracy requirement is undemanding, second-order finite differences may be faster. For the specific case of Poisson's equation on a rectangle, there are specialized fast solution methods that are beyond the scope of this discussion.
## Exercises
1. ✍ Using general $2\times 2$ matrices, verify identities 4 and 6 in {numref}`Theorem {number} `.
(problem-laplace-symmetry)=
2. ✍ Prove that the matrix $\mathbf{A}$ defined in {eq}`poiskronpde` is symmetric if $\mathbf{D}_{xx}$ and $\mathbf{D}_{yy}$ are symmetric.
(problem-laplace-usage)=
3. ⌨ Use {numref}`Function {number} ` to solve the following problems on $[0,1]\times[0,1]$ using $m=n=50$. In each case, make a plot of the solution and a plot of the error.
**(a)** $u_{xx}+u_{yy} = 2x^2[x^2(x-1)+(10x-6)(y^2-y)]$, with $u=0$ on the boundary.
Solution: $u(x,y) = x^4(1-x)y(1-y)$.
**(b)** $u_{xx}+u_{yy} = \left(16 x^2 + (1-4 y)^2\right) \sinh (4 x y-x)$, with $u=\sinh(4xy-x)$ on the boundary.
Solution: $u(x,y) = \sin(4\pi x))$.
**(c)** $u_{xx}+u_{yy} = -(20\pi^2) \sin (4\pi x) \cos (2\pi y)$, with $u = \sin (4\pi x) \cos (2\pi y)$ on the boundary.
Solution: $u(x,y) = \sin (4\pi x) \cos (2\pi y)$.
4. ⌨ For each case in Exercise 3, solve the problem using {numref}`Function {number} ` with $m=n=20,30,40,\ldots,120$. For each numerical solution compute the maximum absolute error on the grid. On a log-log plot, compare the convergence of the error as a function of $n$ to theoretical second-order accuracy.
5. ⌨ Copy {numref}`Function {number} ` to a new function named `poischeb`, and modify it to use a Chebyshev discretization rather than finite differences. For each item in [Exercise 3](problem-laplace-usage), solve the problem using `poischeb` for $m=n=10,15,20,\ldots,40$. For each numerical solution compute the maximum absolute error on the grid. Show the convergence of the error as a function of $n$ on a log-linear plot.
6. ⌨ Sometimes boundary conditions are specified using a piecewise definition, with a different formula for each side of the domain. Use {numref}`Function {number} ` with $m=n=60$ to solve the Laplace equation on $[0,1]^2$ with boundary conditions
$$
u(0,y) = u(1,y) \equiv 0, \quad u(x,0) = \sin(3\pi x), \quad u(x,1) = e^{2x}(x-x^2).
$$
Make a surface plot of your numerical solution.