ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Mon, 07 Mar 2011 16:30:41 +0100Porting a finite-differences-matrix from Matlab to Numpyhttps://ask.sagemath.org/question/7982/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!Thu, 03 Mar 2011 21:23:00 +0100https://ask.sagemath.org/question/7982/porting-a-finite-differences-matrix-from-matlab-to-numpy/Answer by Jason Grout for <p>Hi there,</p>
<p>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).</p>
<p>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 :)</p>
<p>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).</p>
<p>For starters, here's the Matlab initialization code (my comments):</p>
<pre><code>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
</code></pre>
<p>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.</p>
<pre><code>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;
</code></pre>
<p>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.</p>
<p>If I didn't make myself clear, please ask in the comments.</p>
<p>Thanks!</p>
https://ask.sagemath.org/question/7982/porting-a-finite-differences-matrix-from-matlab-to-numpy/?answer=12176#post-id-12176This page might be useful: [http://www.scipy.org/NumPy_for_Matlab_Users](http://www.scipy.org/NumPy_for_Matlab_Users)Mon, 07 Mar 2011 16:30:41 +0100https://ask.sagemath.org/question/7982/porting-a-finite-differences-matrix-from-matlab-to-numpy/?answer=12176#post-id-12176Answer by cswiercz for <p>Hi there,</p>
<p>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).</p>
<p>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 :)</p>
<p>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).</p>
<p>For starters, here's the Matlab initialization code (my comments):</p>
<pre><code>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
</code></pre>
<p>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.</p>
<pre><code>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;
</code></pre>
<p>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.</p>
<p>If I didn't make myself clear, please ask in the comments.</p>
<p>Thanks!</p>
https://ask.sagemath.org/question/7982/porting-a-finite-differences-matrix-from-matlab-to-numpy/?answer=12169#post-id-12169It 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$ matrix`A` 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](http://docs.scipy.org/doc/numpy/user/basics.html) documentation for an introduction to using these data types.Fri, 04 Mar 2011 13:35:31 +0100https://ask.sagemath.org/question/7982/porting-a-finite-differences-matrix-from-matlab-to-numpy/?answer=12169#post-id-12169Comment by sebastian_k for <p>It seems like your question has three parts: </p>
<ol>
<li>What does this Matlab code do?</li>
<li>How do I turn this into Sage/Numpy/Scipy?</li>
<li>What is this Matlab code used for?</li>
</ol>
<p>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.</p>
<p>The Matlab function <code>mACfd</code> creates an $N_2^2 \times N_2^2$ matrix<code>A</code> and a vector <code>C</code> of $N_2^2$ components. <code>A</code> 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$</p>
<p>Python/Sage/Numpy/Scipy does loops, of course, so starting with the <code>numpy.array</code> or <code>numpy.matrix</code> object you should be on your way to writing the same routine. See the <a href="http://docs.scipy.org/doc/numpy/user/basics.html">Numpy Basics</a> documentation for an introduction to using these data types.</p>
https://ask.sagemath.org/question/7982/porting-a-finite-differences-matrix-from-matlab-to-numpy/?comment=22021#post-id-22021Alright, 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!Fri, 04 Mar 2011 17:37:54 +0100https://ask.sagemath.org/question/7982/porting-a-finite-differences-matrix-from-matlab-to-numpy/?comment=22021#post-id-22021Comment by sebastian_k for <p>It seems like your question has three parts: </p>
<ol>
<li>What does this Matlab code do?</li>
<li>How do I turn this into Sage/Numpy/Scipy?</li>
<li>What is this Matlab code used for?</li>
</ol>
<p>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.</p>
<p>The Matlab function <code>mACfd</code> creates an $N_2^2 \times N_2^2$ matrix<code>A</code> and a vector <code>C</code> of $N_2^2$ components. <code>A</code> 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$</p>
<p>Python/Sage/Numpy/Scipy does loops, of course, so starting with the <code>numpy.array</code> or <code>numpy.matrix</code> object you should be on your way to writing the same routine. See the <a href="http://docs.scipy.org/doc/numpy/user/basics.html">Numpy Basics</a> documentation for an introduction to using these data types.</p>
https://ask.sagemath.org/question/7982/porting-a-finite-differences-matrix-from-matlab-to-numpy/?comment=22022#post-id-22022Sorry, 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! :)Fri, 04 Mar 2011 16:46:26 +0100https://ask.sagemath.org/question/7982/porting-a-finite-differences-matrix-from-matlab-to-numpy/?comment=22022#post-id-22022Comment by cswiercz for <p>It seems like your question has three parts: </p>
<ol>
<li>What does this Matlab code do?</li>
<li>How do I turn this into Sage/Numpy/Scipy?</li>
<li>What is this Matlab code used for?</li>
</ol>
<p>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.</p>
<p>The Matlab function <code>mACfd</code> creates an $N_2^2 \times N_2^2$ matrix<code>A</code> and a vector <code>C</code> of $N_2^2$ components. <code>A</code> 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$</p>
<p>Python/Sage/Numpy/Scipy does loops, of course, so starting with the <code>numpy.array</code> or <code>numpy.matrix</code> object you should be on your way to writing the same routine. See the <a href="http://docs.scipy.org/doc/numpy/user/basics.html">Numpy Basics</a> documentation for an introduction to using these data types.</p>
https://ask.sagemath.org/question/7982/porting-a-finite-differences-matrix-from-matlab-to-numpy/?comment=22020#post-id-22020That's good to hear!Fri, 04 Mar 2011 17:50:29 +0100https://ask.sagemath.org/question/7982/porting-a-finite-differences-matrix-from-matlab-to-numpy/?comment=22020#post-id-22020Comment by cswiercz for <p>It seems like your question has three parts: </p>
<ol>
<li>What does this Matlab code do?</li>
<li>How do I turn this into Sage/Numpy/Scipy?</li>
<li>What is this Matlab code used for?</li>
</ol>
<p>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.</p>
<p>The Matlab function <code>mACfd</code> creates an $N_2^2 \times N_2^2$ matrix<code>A</code> and a vector <code>C</code> of $N_2^2$ components. <code>A</code> 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$</p>
<p>Python/Sage/Numpy/Scipy does loops, of course, so starting with the <code>numpy.array</code> or <code>numpy.matrix</code> object you should be on your way to writing the same routine. See the <a href="http://docs.scipy.org/doc/numpy/user/basics.html">Numpy Basics</a> documentation for an introduction to using these data types.</p>
https://ask.sagemath.org/question/7982/porting-a-finite-differences-matrix-from-matlab-to-numpy/?comment=22023#post-id-22023I 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 http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.lstsq.html#numpy.linalg.lstsqFri, 04 Mar 2011 16:07:04 +0100https://ask.sagemath.org/question/7982/porting-a-finite-differences-matrix-from-matlab-to-numpy/?comment=22023#post-id-22023Comment by sebastian_k for <p>It seems like your question has three parts: </p>
<ol>
<li>What does this Matlab code do?</li>
<li>How do I turn this into Sage/Numpy/Scipy?</li>
<li>What is this Matlab code used for?</li>
</ol>
<p>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.</p>
<p>The Matlab function <code>mACfd</code> creates an $N_2^2 \times N_2^2$ matrix<code>A</code> and a vector <code>C</code> of $N_2^2$ components. <code>A</code> 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$</p>
<p>Python/Sage/Numpy/Scipy does loops, of course, so starting with the <code>numpy.array</code> or <code>numpy.matrix</code> object you should be on your way to writing the same routine. See the <a href="http://docs.scipy.org/doc/numpy/user/basics.html">Numpy Basics</a> documentation for an introduction to using these data types.</p>
https://ask.sagemath.org/question/7982/porting-a-finite-differences-matrix-from-matlab-to-numpy/?comment=22024#post-id-22024Thanks 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?Fri, 04 Mar 2011 15:51:28 +0100https://ask.sagemath.org/question/7982/porting-a-finite-differences-matrix-from-matlab-to-numpy/?comment=22024#post-id-22024