Developing ideas in Matlab : You are doing it right!

Hi all!

I have collected a bunch of things that I want to write about. The graduate courses I am taking are full of interesting stuff, some of which I haven’t seen before and I will start commenting on them the following days. Today I want to say “thank you” to Mathworks for providing us with a tool on which we can develop ideas in a much faster and easier way, keeping our mind more over the details and not over the implementation of the code.

But let’s be specific. One of the first homeworks that I had to give was about finite differences and fast Poisson solvers (which is a pretty cool concept, I am gonna write something about it in the following days). The idea was that we had a Poisson equation in [0,1], u_{xx}=f, u(1)=u(0)=1, with f being a somewhat complex function, and we had to solve this using finite differences. Let’s write it in a slightly more general way.

u_{xx}=F,\ u(1)=a,\ u(0)=b,\ a,b\in\mathbb{R}

Now I had coded finite differences again in the past. In fact, I had done it in Matlab but because the professor which taught the course didn’t want me to use any function already written in Matlab, I had to code everything myself.

So, what I had to do was code the construction of a sparse, banded matrix A, write the code which would do the Cholesky decomposition of the matrix and then use backward and forward substitution to solve the system. Oh and I forgot to mention that it was my code that did the forward and backward substitution. So, while pedagogical, the whole thing took an awful time to write and any errors that happened (which are bound to happen  when writing many lines of code) further distracted me from the point, which was to solve the Poisson equation. All in all, I did it in around 120 lines of code.

You can imagine my happiness when my now numerical solution of PDEs professor told me that I can use whatever I wanted. Ladies and gentlemen, here is a way to implement finite differences in Matlab, using m+2 points in [0,1] :

x = ((1 : (m+2)) - 1)*1/(m+1);
F = f(x);
A = (m+1)^2*sparse([1:m,2:m,1:m-1],[1:m,1:m-1,2:m],[-2*ones(1,m),1*ones(1,m-1),1*ones(1,m-1)]);
F(1) = F(1) - (m+1)^2*a;
F(end) = F(end) - (m+1)^2*b;
Uappr = [a,(A\F.').',b];

And it only took six lines of code! Well, somebody might want to add a little more to make it a function and perhaps in order to produce graphs or compute errors, but this code does exactly what it should! Feed it with a right hand side function f(x) and a number of nodes m and bam! You ‘ve got your (approximate) solution. I know that the nodes seem m+2 in the code but there isn’t any reason to use the first and last nodes. We have good boundary conditions and already know the solution at those points.

So, I was able to write and test my code fast and I know that it works OK. Since my first concern wasn’t speed but to understand what I was doing, using Matlab was a good idea. If I want something more specialized or hand-tuned, I can convert it to C or Fortran and optimize it further. That goes to all die-hard “I won’t use Matlab cause it’s not a programming language” fans. It’s easier to find all the logical errors of your code if you don’t have to worry that much about all the extra stuff that come with lower lever languages, i.e. memory allocation, use of an efficient data structure and finding the right library for the job. 🙂

Until next time.