I need to examine the (static) electric field distribution around various electrode configurations in the presence of a dielectric, and have been using a finite difference approach to the 2D generalised Poisson equation do this. James Nagel's nice paper at http://www.drjamesnagel.com/notes/Nagel%20-%20Numerical%20Poisson.pdf has been my guide and suggests SOR on a discretized integral form of Poisson's equation. But Nagel's paper also notes that for large dielectric constants - which I have - SOR converges slowly or fails: then he uses the Matlab/Octave "\" as a blackbox operation on the sparse system matrix. This works fine for me on problems up to 300x300 grid or thereabouts - but it is slow and I really want to use C or something similar to integrate with other routines. So I am seeking advice on a suitable algorithm to apply to this problem. I think ADI may be what I need - but before learning how to code and use it, I'd value guidance and perhaps alternative recommendations. Perhaps I just need to call the LAPACK version of "\"? Of course, I want unconditional stability, high accuracy and instant convergence! Thanks for reading!
There are many excellent codes out there that solve the Poisson equation. I may point you at step-6 of the deal.II tutorial at http://dealii.org/developer/doxygen/deal.II/step_6.html as just one example. (Disclaimer: That's my own code.) Most other finite element libraries will have similar example codes.
A good way to approach this is by using existing powerful numerical packages such as PETSc. Look for PETSc examples of solving the Poisson equation, and that will be a good starting point.
The FDM paper (which I personally wrote, BTW) is intended more for basic education than it is for large-scale implementations. As you may have noticed, the major bottleneck for large domains is setting up the system matrix and inverting it. In my experience, Matlab simply runs out of memory when trying to run the inversion outright, so you have to be clever how this step is implemented. I did my best to try and streamline this step in the code, but getting a solid matrix inversion for this kind of computation is a numerical challenge unto itself. If anyone has some clever updates for allowing my codes to process larger matrices, I'd love to hear about it.
... or of course you can always just use FEM instead, and simply sample the grid more efficiently. :)
EDIT: I'm not sure about the differences between Matlab and Octave when dealing with sparse matrices. In Matlab, setting up the matrix has always been the easy part for me. It has always been the inversion step that kills it by running out of memory.
I suggest using the "Run and Time" option in MATLAB to see exactly where the time is being spent in your 300x300 grid example. I am almost certain it is not in the \ operator. Generally, MATLAB has very fast, state-of-the-art implementations for sparse matrix factorizations.
I took at look at James Nagel's MATLAB code at the site referred to in your link. I didn't try to run it since I couldn't find an example case. But I did notice that there is a loop over all points in the grid (e.g. 300x300) to fill the sparse matrix prior to calling . I strongly suspect this is where most of the time is being spent. It is generally more efficient in MATLAB to write the code so that operations are performed on long vectors rather than scalars.
You can find many references on this topic if you search for "vectorized MATLAB code".
I should add that if you are interested in a fast, basic finite element code for the Poisson equation written in MATLAB, you can find one at the following link
along with a discussion of how it was designed.
So, finding a good solution has been instructive. First point is that as Bill suggested, setting up the system in Octave (rather than solving it) takes most of the time - but it is not easy code to vectorise. Also, for my problem formulation with mixed boundary conditions, the system matrix is not symmetric, which rules out some efficient approaches. I looked at different libraries and now use the sparse LU solver in umfpack, included in my Linux system as part of the suitesparse library. I gather both Octave and Matlab use umfpack which I take as a recommendation - it is a traditional library, easy to call from C as well as C++, and has excellent, clear documentation. A straightforward C routine solves a 320x320 FDE grid >100 times faster than in Octave and easily copes with 1000x1000 in 30s on a cheap laptop.
FYI, I've recently published a more comprehensive tutorial with the IEEE APS magazine. Here is the direct link. This one includes more material and has been formally peer-reviewed, so publication references finally have something solid to point to.
It's a great starter paper to get the basic idea of FDM working. It also seems to work very well with domain sizes on the order of 300x300 or so, depending on your memory. As we've all seen, the only real bottleneck for larger domains is setting up that pesky system matrix and inverting it.