As part of an assignment, I had to write Matlab code to solve Poisson’s equation in a trapezoidal domain, by transforming the equation and the BCs to the unit square. The transformed PDEs were ugly, and figuring out the structure of the associated system matrix didn’t seem like fun, so I wrote up some pretty neat SICP-inspired utility code that makes coding finite differences (on a uniform grid in a square domain) pretty darn easy.
The idea behind it is that the real devil in writing fd code is that you have transform a grid function (say you’re solving for 
on a square domain) into a vector so you can use linear algebra to solve for it; your neighbor relations go all wacky when you do this, so it’s hard to write out the system matrix that arises. The best way, IMHO, to defeat this devil is to delegate responsibility for coordinate management– write everything in terms of the original grid coordinates and let the computer handle the translation to vector coordinates (I’ve heard they’re good at those sort of brainless tasks).
At the heart of the system are two functions:
-
indices = gridtovect(X, Y, N)– which given vectors of grid x-coordinates X = [x1, x2, ...] and grid y-coordinates Y= [y1, y2, ...] and a gridsize N returns the linear indices of the points (x1, y1), (x2, y2), ... (using the transformation 
).
-
indices = tensorind(X, Y, N)– which given vectors of system matrix coordinates X = [x1, x2, ...] and Y= [y1, y2, ...] and a gridsize N, returns the linear indices of the system matrix coordinates (x1, y1), (x2, y2), ... (recall in Matlab you can treat a matrix of arbitrary dimension as a linear vector; the index transformations are done with ind2sub and sub2ind)
For convenience these functions will assume if you pass in a constant in place of a vector, that you mean a constant vector of appropriate size. Using these, I constructed functions which modify the system matrix, adding a system of equations that all of the indicated points of the grid must satisfy, e.g.:
[code]
% given a system matrix A, and x,y pairs in column vectors X,Y, modifies
% the entries in A to add a 2nd order accurate centered finite difference
% term to approximate the x-deriv, centered at each x,y pair; accepts a
% vector field c to multiply each FD by, and the step size dx
function A = centerFDd1x(A, X, Y, dx, c, N)
self = gridtovect(X, Y, N);
lft = gridtovect(X-1, Y, N);
rt = gridtovect(X+1, Y, N);
A( tensorind(self, lft, N) ) = …
A( tensorind(self, lft, N) ) + c * -1/(2*dx);
A( tensorind(self, rt, N) ) = …
A( tensorind(self, rt, N) ) + c * 1/(2*dx);
end
[/code]
Using my problem as an example, the pde that had to be satisfied by the points on the interior of the square was
Here’s how I added the entries for the interior nodes to the system matrix (the names of x and y were changed to xi and eta to protect the innocent to serve as a reminder that this is in a transformed domain):
[code]
%% setup the interior PDE
%% u_{xi, xi} c1(xi, eta) + u_{xi, eta} c2(xi, eta)
%% + u_{xi} c3(xi,eta) + u_{eta, eta} = -h^2
xi = dxi:dxi:1-dxi;
for grideta = 2:N-1
eta = (grideta -1)*dxi;
c1 = 4*(h^2 + a^2*xi.^2) / (b + 2*a*eta)^2;
c2 = -(4*a*xi) / (b + 2*a*eta);
c3 = 8*a^2*xi / (b + 2*a*eta)^2;
% u_{xi, xi} c1(xi, eta)
A = centerFDd2x(A, 2:N-1, grideta, dxi, c1, N);
% u_{xi, eta} c2(xi, eta)
A = asymFDdxy(A, 2:N-1, grideta, dxi, c2, N);
% u_{xi} c3(xi, eta)
A = centerFDd1x(A, 2:N-1, grideta, dxi, c3, N);
% u_{eta, eta}
A = centerFDd2y(A, 2:N-1, grideta, dxi, 1, N);
self = gridtovect(2:N-1, grideta, N);
f(self) = -h^2;
end
[/code]
Check out the code if you’re interested (the top stuff was for my particular problem, the reusable code’s at the bottom).
