Ask Your Question

Porting a finite-differences-matrix from Matlab to Numpy

asked 2011-03-03 14:23:00 -0500

sebastian_k gravatar image

Hi there,

I'm an engineering student plagued by a prof who adores Matlab. For reasons I'm sure you understand, I'd like to avoid Matlab and learn Numpy/Scipy instead. Now I'm given this code, which I'm supposed to adapt to a problem given in class (irrelevant).

If you're somewhat familiar with Matlab and interested in helping me port a weird piece of code, then this question is for you -- read on :)

Short explanation: We're creating a meshgrid of nodes, and assign voltages (matrix Vstart) to some of them. We'd like to find the potentials of the nodes in between. You'd assume we want the voltage drop to follow ?²=0, and indeed I get a fine, smooth potential drop if I loop through the in-between-nodes 50 times and repeatedly assign to them the average value of the nodes around them (that was my approach -- code not shown here).

For starters, here's the Matlab initialization code (my comments):

for m = 1:length(N)
d = a/N(m);                  #distance between mesh nodes
N1 = N(m)+1;                 #side length of the central conductor
N2 = b/a*N(m)+1;             #side length of meshgrid

Vstart = zeros(N2,N2);
Vstart(1,:) = Vb;            #set the outer nodes to Vb
Vstart(:,1) = Vb;
Vstart(N2,:) = Vb;
Vstart(:,N2) = Vb;
Vstart(lim1:lim2,lim1:lim2) = Va;    #set the inner nodes to Va

Now we're told to use a different approach than the one described above, and we're given a piece of Matlab code that fills a matrix A with a bunch of negative-quarters, and another matrix A with a long line of 0's, Vb's and Va's. It then solves a system and puts the results into V. I'm not quite sure how the matrix works here, although the many 1/4's make me think this is just my previous approach in disguise.

function[A,C] = mACfd(Vstart,N2)
C = zeros(N2^ 2,1);
for i = 1:N2
for j = 1:N2
k = (i-1)*N2+j;
A(k,k) = 1;
if (Vstart(i,j) == 0)
A(k,k-N2) = -1/4;
A(k,k+N2) = -1/4;
A(k,k+1) = -1/4;
A(k,k-1) = -1/4;
else C(k) = Vstart(i,j);

[A,C] = mACfd(Vstart,N2);
V = inv(A)*C;
for i = 1:N2
V2D(i,:) = V((i-1)*N2+1:i*N2);

I don't understand how matrices A and C are filled, what their dimensions are, what they're for, and how I can efficiently reproduce the algorithm in Python. Particularly the weird indexing irks me.

If I didn't make myself clear, please ask in the comments.


edit retag flag offensive close merge delete

2 answers

Sort by » oldest newest most voted

answered 2011-03-04 06:35:31 -0500

updated 2011-03-07 09:12:26 -0500

It seems like your question has three parts:

  1. What does this Matlab code do?
  2. How do I turn this into Sage/Numpy/Scipy?
  3. What is this Matlab code used for?

This first question is perhaps not best asked in this forum but I'll provide some insight. For the second question I'll suggest some outside resources that can get you started. As for that last question, I'm afraid you'll have to ask your professor about that.

The Matlab function mACfd creates an $N_2^2 \times N_2^2$ matrixA and a vector C of $N_2^2$ components. A is a banded matrix. (i.e. Tridiagonal with additional diagonals distance $N_2$ away.) The way Matlab is creating this matrix-vector pair is by means of a nested for loop across the indices $i$ and $j$ with an auxillary index $k = (i-1)N_2+j$

Python/Sage/Numpy/Scipy does loops, of course, so starting with the numpy.array or numpy.matrix object you should be on your way to writing the same routine. See the Numpy Basics documentation for an introduction to using these data types.

edit flag offensive delete link more


Thanks Chris, I understand how it works now -- how the matrix is created, how it works, and how I can do the same thing in sage/numpy. I have a problem inverting/solving/lstsq'ing though: as soon as the number of equations exceeds 160, sage segfaults. I tried this example:, and again, setting Nparam to anything higher than 160 causes a segfault. (numpy 1.3.0). I'm now compiling 1.5.0 and hope that'll solve it; if not, any idea what could cause this error?

sebastian_k gravatar imagesebastian_k ( 2011-03-04 08:51:28 -0500 )edit

I executed the example code you mentioned (the lstsq example) and did not get any segfaults. Which version of Sage are you running? I rand the code on version 4.6.1 which comes with Numpy version 1.5.0. (Run "import numpy; print numpy.__version__" in Sage.) Finally, I'm not familiar with how to use Numpy's lstsq function. You may want to read its documentation very closely at

cswiercz gravatar imagecswiercz ( 2011-03-04 09:07:04 -0500 )edit

Sorry, the sample code works fine as it is -- only if I increase Nparam to, say, 200 I get a segfault (within ipython, that is; the notebook just doesn't give me output). That was using numpy 1.3.0. I am now upgrading sage to numpy 1.5.0 and don't get the error anymore, although it's still not working as I'd like it to. When it's done compiling all the upgrades I'll experiment some more and let you know how it goes. A thousand thanks for your help! :)

sebastian_k gravatar imagesebastian_k ( 2011-03-04 09:46:26 -0500 )edit

Alright, it works fine now with numpy 1.5.0, the bug disappeared. Problem solved. Now I can start pythonifying -- the code is already shorter than its Matlab equivalent :) Thanks again!

sebastian_k gravatar imagesebastian_k ( 2011-03-04 10:37:54 -0500 )edit

That's good to hear!

cswiercz gravatar imagecswiercz ( 2011-03-04 10:50:29 -0500 )edit

answered 2011-03-07 09:30:41 -0500

Jason Grout gravatar image
edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools


Asked: 2011-03-03 14:23:00 -0500

Seen: 893 times

Last updated: Mar 07 '11