Firstly, I'm working with Dirichlet boundary conditions, and I arbitrarily set the edge of my matrix to have these constant boundary conditions throughout the problem. There is a heat-generating source in the middle of the world, and there is some geometry throughout the world, and of course, at the edge of the matrix (that is, at the edge of the world), the temperature is always 293K. Air is the conducting medium between the geometry and the edge of the universe.
I noticed that the del2.m file is broken up into three sections:
- initialization and option management
- calculation of the central areas of the matrix
- calculation of the edges
Since the edges aren't necessary (that is, they represent the boundary conditions), I ignore the last section. Also, I'm working only with 2-d matrices right now.
So to calculate the central areas of the matrix, I need an array that computes the del2 and another which sets the edges to zero:
## Declare the derivative matrix
dmat = diag(-1*ones(1,Nux-1),1) + 2*eye(Nux) + \
diag(-1*ones(1,Nux-1),-1);
dmat = dmat ./ -4.0;
## Declare the blanking matrix
blankmat = blkdiag([0],ones(Nux-2,Nuy-2),[0]);
Next, here's the calculation for my del2:
( (x*dmat)+(dmat*x) ) .* blankmat / dx / dy
So I go ahead and run this in my heat model and get a picture that looks exactly like what I was getting with the del2 operator after 50 cycles. I also calculated the difference, and only the area adjacent to the edges is different, and that only by 0.001K, which is quite acceptable and certainly within numerical error (representing an error of 3.4e-06).
The increase in speed is about 25x over using del2.m