The goal is to solve\[ \frac {\partial ^{2}u}{\partial x^{2}}+\frac {\partial ^{2}u}{\partial y^{2}}+\frac {\partial ^{2}u}{\partial z^{2}}=-f(x,y,z) \] On the unit cube. The following diagram was made to help setting up the 3D scheme to approximate the above PDE
The discrete approximation to the Laplacian in 3D is\[ \frac {\partial ^{2}u}{\partial x^{2}}+\frac {\partial ^{2}u}{\partial y^{2}}+\frac {\partial ^{2}u}{\partial z^{2}}=\frac {1}{h^{2}}\left (U_{i-1,j,k}+U_{i+1,j,k}+U_{i,j-1,k}+U_{i,j+1,k}+U_{i,k,k-1}+U_{i,j,k+1}-6U_{i,j,k}\right ) \] For the direct solver, the \(A\) matrix needs to be formulated. From\[ \frac {1}{h^{2}}\left ( U_{i-1,j,k}+U_{i+1,j,k}+U_{i,j-1,k}+U_{i,j+1,k}+U_{i,k,k-1}+U_{i,j,k+1}-6 U_{i,j,k}\right ) =f_{i,j,k}\] Solving for \(U_{i,j,k}\) results in\[ U_{i,j,k}=\frac {1}{6}\left ( U_{i-1,j,k}+U_{i+1,j,k}+U_{i,j-1,k}+U_{i,j+1,k}+U_{i,k,k-1}+U_{i,j,k+1}-h^{2}f_{i,j,k}\right ) \] To help make the \(A\) matrix, a small example with \(n=2,\) is made. The following diagram uses the standard numbering on each node
By traversing the grid, left to right, then inwards into the paper, then upwards, the following \(A\) matrix results
The recursive pattern involved in these \(A\) matrices can now be seen. Each \(A\) matrix contains inside it a block on its diagonal which repeats \(n\) times. Each block in turn contain inside it, on its diagonal, smaller block, which also repeats \(n\) times.
It is easier to see the pattern of building \(A\) by using numbers for the grid points, and label them in the same order as they would be visited, this allowes seeing the connection between each grid point to the other easier. For example, for \(n=2\),
The connections are now more easily seen. Grid point 1 has connection to only points \(2,3,5\). This means when looking at the \(A\) matrix, there will be a \(1\) in the first row, at columns \(2,3,5\). Similarly, point \(2\) has connections only to \(1,4,6\), which means in the second row there will be a \(1\) at columns \(1,4,6\). Extending the number of points to \(n=3\) to better see the pattern of \(A\) results in
From the above it is seen that for example point \(1\) is connected only to \(2,4,10\) and point \(2\) is connected to \(1,3,5,11\) and so on.
The above shows that each point will have a connection to a point which is numbered \(n^{2}\) higher than the grid point itself. \(n^{2}\) is the size of the grid at each surface. Hence, the general \(A\) matrix, for the above example, can now be written as
The recursive structure can be seen. There are \(n=3\) main repeating blocks on the diagonal, and each one of them in turn has \(n=3\) repeating blocks on its own diagonal. Here \(n=3\), the number of grid points along one dimension.
Now that the \(A\) structure is understood, the A matrix can be generated.
Example 1: Using \(n_x=2\), \(n_y=2\), \(n_z=2\). These are the number of grid points in the \(x,y,z\) directions respectively.
Example 2: Using \(n_x=2\), \(n_y=2\), \(n_z=3\).
Example 3: Using \(n_x=3\), \(n_y=3\), \(n_z=3\).
Matlab
nx=3; ny=3; ex=ones(nx,1); Lx=spdiags([ex -3*ex ex],[-1 0 1],nx,nx); ey=ones(ny,1); Ly=spdiags([ey -3*ey ey],[-1 0 1],ny,ny); Ix=speye(nx); Iy=speye(ny); L2=kron(Iy,Lx)+kron(Ly,Ix); nz=3; N=nx*ny*nz; e=ones(N,1); L=spdiags([e e],[-nx*ny nx*ny],N,N); Iz=speye(nz); A=kron(Iz,L2)+L; full(A)
