I want to solve the 2D Laplace's equation: $$ \frac{\partial^2 T}{\partial x^2 } + \frac{\partial^2 T}{\partial y^2 } = 0 $$ with boundary conditions: T(x=0)=T(x=1)=T(y=1)=0 and T(y=0)=1 on a 1024x1024 discrete grid.
How do I do this?
Is it the Discrete Sine Transform I should use? If so which of them?
So many options!
This can be solved analytically using separation of variables. So that's good. Gives us something to check against, make sure that we are doing this right, and a good way of checking convergence.
The non-continuous aspect isn't a problem. Might need a few more terms to converge, but will still converge.
1024x1024 I think will be overkill, to be honest, depending of course on your required accuracy, even a quick and dirty finite difference solution will give good results there. No big deal, though whatever you do, I highly recommend coding in general and setting the grid to be much coarser while testing. Like 2x2 or 3x3. That makes it much easier to see when (probably not if) you mess something up and make it explode, not to mention that it's much faster for the computer to solve.
Ok, enough preamble, basic idea: you choose a series solution to approximate the true solution and determine the coefficients that best satisfy the problem. Let's look how I would approach this problem. First, let's think about what we want our solution to look like:
$\Theta(x,y)=\sum_{i=0}^n\sum_{j=0}^ma_{i,j}F_i(x)G_j(y)$
(Note, I'm using $\Theta(x,y)$ to denote our approximate solution to the true solution $T(x,y)$) We have to choose the $F_i$'s & $G_j$'s (our basis functions) and determine the $a_{i,j}$'s that best solve the problem.
Now, you can choose anything you want for $F_i$'s & $G_j$'s, as long as they are all linearly independent (orthogonality is nice, but not strictly required). So let's choose basis functions that make our life easier. For a problem like this, I would chose basis functions that automatically satisfy the boundary conditions. So in our $x$-coordinate, let's use:
$F_i(x) = \sin(i\pi x)$ for $i=1...n-1$
which is $0$ at at $x=0$ and $x=1$ for all $i$, and independent of $G_j$, so those BC's are taken care of for $x$.
The $y$-coordinate is a little trickier (if we use $\sin (j\pi x)$, we can't satisfy our non-homogeneous BC), but not too bad:
$G_j(y) = \cos((j+1/2)\pi y)$ for $j=1...m$
This satisfies our BC of $T=0$ at $y=1$. And $G_j(0) = 1$ for all $j$, so that's something we can work with.
Let's look back at our desired equation:
$\Theta(x,y)=\sum_{i=0}^n\sum_{j=0}^ma_{i,j}F_i(x)G_j(y) = \sum_{i=1}^{n-1}\sum_{j=1}^ma_{i,j}\sin(i\pi x)\cos((j+1/2)\pi y)$
(Note: since we've handled the 3 of the 4 homogeneous BC, I've adjusted the limits of the summation to go from $1$ to $n-1$ and $1$ to $m$ rather than $0$ to $n$ and $0$ to $m$ )
Okay, now we need to find the coefficients, $a_{i,j}$, which give the best solution. Well, first, we need to define what we mean by "best". Ideally, we want our approximate solution $\Theta(x,y)$ to equal the real solution $T(x,y)$, or at least minimize $\Theta(x,y)-T(x,y)$ in some averaged sense. But we don't know what $T(x,y)$ is, in general, all we have to work with is the governing equation (Laplace's Eq, in this case). In other words, we will chose our $a_{i,j}$ in such a way that $\Theta(x,y)$ approximately satisfy the governing equation, so let's plug our solution into that
$\frac{\partial^2 \Theta}{\partial x^2} + \frac{\partial^2 \Theta}{\partial y^2} = 0$
$\sum_{i=1}^{n-1}\sum_{j=1}^m-i^2\pi^2 a_{i,j}\sin(i\pi x)\cos((j+1/2)\pi y) - \sum_{i=1}^{n-1}\sum_{j=1}^m(j+1/2)^2\pi^2a_{i,j}\sin(i\pi x)\cos((j+1/2)\pi y) = 0$
In general, we cannot choose a set of $a_{i,j}$ that ensures that this equation is satisfied everywhere (i.e. for all $x$ and $y$), and the deviation from the solution is called the "residual". The "best" solution would be to use the Galerkin approach to minimize the sum of square of the residual across the entire domain. But it involves performing integrations and can get rather messy and involved in the general case (though for this case, it may be easier, since you can take advantage of the orthogonality of trig functions and everything is nice and linear).
Alternatively, you can just set the residual to be zero at $(n-1)$ in discrete points $x$ and $(m-1)$ in discrete points $y$. This is called collocation tends to be much simpler than the Galerkin approach. You can use equal spaced nodes using this technique, but you'll get better results if you use Chebyshev nodes (that will increase node density near the boundaries, and prevent Runge's phenomena from giving you trouble).
Either approach will give you $(n-1)\times (m-1)$ linear equations to solve simultaneously for the $((n-1)\times m)$ $a_{i,j}$'s. Obviously, that's not enough equations, but remember, we still have that final, non-homogeneous BC to take care of:
$\Theta(x,0)= \sum_{i=1}^{n-1}\sum_{j=1}^ma_{i,j}\sin(i\pi x)=1$
We can use this to provide the remaining necessary constraints to fully solve for our coefficients. Again, we can use collocation to specify $\Theta(x,0)=1$ at $(n-1)$ specific values of $x$ or we can use a Galerkin approach (which would be a similar procedure to a Fourier transform in this case). Either way, you'll get another $(n-1)$ equations to work with, so you have $(n-1)\times (m)$ linear equations in $(n-1)\times m$ unknowns. And it's just a matter of using your linear solver of choice.
There's also a technique called the "Boundary Element Method", which I think is really cool and something you may be interested in. Basically you only discretize and solve along the boundary and can calculate the solution at any point on the interior after the fact. It effectively reduces the dimensionality of the problem by one (i.e. a $2D$ becomes $1D$) and is especially useful where the ratio of the boundary area to the volume is small. It's limited in application though, I think the governing equation has to be homogeneous, for example, but it's very powerful when it can be used.
