sebastian_k's profile - activity

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.

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!

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! :)

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: http://www.scipy.org/Numpy_Example_List#head-ef374c752c3530b09022336de9adc2ce85af822c, 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?

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;
lim1=(N2-N1)/2+1;
lim2=(N2+N1)/2;
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);
end;
end;
end;

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

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.

Thanks!