Ask Your Question

Revision history [back]

click to hide/show revision 1
initial version

Porting a finite-differences-matrix from Matlab to Numpy

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!