1

I'm a C and MATLAB user. When I started learning Python (a week ago) I noticed that I don't use full potential of MATLAB, in particular array operations. I use for loops often, probably because I learnt programming in C.

In a previous tip, I learnt to use cumsum and other efficient array operations, for example:

alpha = [1e-4,1e-3,1e-4,1e-1,1e-2,1e-3,1e-6,1e-3];
zeta = alpha / (dz*dz)
nz = 101
l=[0.3,0.1,0.2,0.1,0.1,0.1,0.2];
wz = cumsum(l*(nz-1));
nl = lenght(l);

Is it possible simplify the following code in Python (Numpy) or MATLAB?

      A = zeros(nz,nz);
      i=1;
      for j = 2:wz(i)-1
        A(j,j-1) = zeta(1,1);
        A(j,j) = -2*zeta(1,1);
        A(j,j+1) = zeta(1,1); % layer 1 nodes 
      end

      %cicle to n-layers
      for i=2:nl
          for j=wz(i-1):wz(i-1)
              A(j,j-1) = zeta(1,i-1);
              A(j,j) = -zeta(1,i-1)-zeta(1,i);
              A(j,j+1) = zeta(1,i); 
          end

          for j=wz(i-1)+1:wz(i)
              A(j,j-1) = zeta(1,i);
              A(j,j) = -2*zeta(1,i);
              A(j,j+1) = zeta(1,i);
          end

      end
end
Improve this question
4
  • 2
    It would be easier for us if you would explain in a sentence or two what the code is supposed to do, not just give us the plain code... Commented Jan 16, 2012 at 12:21
  • @Hans This code belongs to an 1D Equation Heat Solver applied to a multilayer. alpha is an array with the diffusivities per layer, l is an array with the heights per layer, wz is an array which aggregates the cumulative sum of points (discrete points) and A is the "state matrix". After calculate state matrix, I'll implement the ode solver. Commented Jan 16, 2012 at 13:10
  • @marco: in case you didn't realise, MATLAB has various ODE solvers built in, while for Python, you can find ODE solvers in Scipy. Commented Jan 16, 2012 at 15:03
  • I'll use a built in ODE solver. But matrix A is an Input of solver. THis code is corret (solves the problem) but it is not efficient. Commented Jan 16, 2012 at 15:12

2 Answers 2

Reset to default
1

I've modified the code below after having a chance to run it on my machine side by side yours. There are still a couple of questions (is A suppose to get larger in the final loop?, what is dz?). The problem you ran into before running this was that I forgot idx_matrix had to be logical.

dz=0.1;
alpha = [1e-4,1e-3,1e-4,1e-1,1e-2,1e-3,1e-6,1e-3];
zeta = alpha / (dz*dz);
nz = 101;
l=[0.3,0.1,0.2,0.1,0.1,0.1,0.2];
wz = cumsum(l*(nz-1));
nl = length(l);

A = zeros(nz);
i=1;

%replaces 1st loop
j_start = 2;
j_end = wz(i)-1;

idx_matrix = false(size(A));
idx_matrix(j_start:j_end,j_start:j_end) = eye(j_end-j_start+1);
A(idx_matrix) = -2*zeta(1,1);

idx_matrix(idx_matrix) = false;
idx_matrix(j_start:j_end,j_start-1:j_end-1) = eye(j_end-j_start+1);
A(idx_matrix) = zeta(1,1);

idx_matrix(idx_matrix) = false;
idx_matrix(j_start:j_end,j_start+1:j_end+1) = eye(j_end-j_start+1);
A(idx_matrix) = zeta(1,1);

%cicle to n-layers
for i=2:nl

    %replaces 3rd loop
    j_start = wz(i-1);
    A(j_start,j_start) = -zeta(1,i-1)-zeta(1,i);
    A(j_start,j_start-1) = zeta(1,i-1);
    A(j_start,j_start+1) = zeta(1,i);

    %replaces 4th loop
    j_start = wz(i-1)+1;
    j_end = min(wz(i),size(A,2)-1);
    idx_matrix = false(size(A));
    idx_matrix(j_start:j_end,j_start:j_end) = eye(j_end-j_start+1);
    A(idx_matrix) = -2*zeta(1,i);

    idx_matrix(idx_matrix) = false;
    idx_matrix(j_start:j_end,j_start-1:j_end-1) = eye(j_end-j_start+1);
    A(idx_matrix) = zeta(1,i);

    idx_matrix(idx_matrix) = false;
    idx_matrix(j_start:j_end,j_start+1:j_end+1) = eye(j_end-j_start+1);
    A(idx_matrix) = zeta(1,i);

end
Improve this answer
Sign up to request clarification or add additional context in comments.

7 Comments

thank you for answering. Your code has an error ??? Subscript indices must either be real positive integers or logicals. A(idx_matrix) = -2*zeta(1,1); It's so early in the code, that I can't figure what you want to do. Yes, third cycle, can be an assignment, I didn't find a better way to this instruction. (This represents the boundary between layers, if you read my code description to @Hans)
AHh, and "you're A matrix" is an array (A=zeros(nz)) is an error or not?
A = zeros(nz) is the same as A = zeros(nz,nz)
Did you intend for A to get larger as this code continues to run? The final for loop adds a column to A because wz(i) == size(A,2) and A(j,j+1) indexes out of bounds, so another column is added. Also, you never defined dz.
you're right. you're code works. However, in the final loop the "zeta" fators are not mulpying per coefficients. Is it possible keep general? If the size of l is chaging?
|
1

To simplify your loops, you can use the function spdiags.

http://www.mathworks.fr/help/techdoc/ref/spdiags.html

For instance your first loop can be written as:

A=full(spdiags(repmat([zeta(1,1),-2*zeta(1,1),zeta(1,1)],wz(i),1),[-1 0 1],wz(i),wz(i)))
Improve this answer

7 Comments

thank you for the reference, but this don't resolve my problem =) I'm trying to find a way to simplify the problem, in other words transform the iteration problem in vector problem
I edited the answer to show you how to use spdiags in your case,
I didn't get it, and when I paste the code in Matlab, returns ??? Index exceeds matrix dimensions. Error in ==> spdiags at 114 a((len(k)+1):len(k+1),:) = [i i+d(k) B(i+(m>=n)*d(k),k)];
works, but the results are different. "My matrix" has a diagonal (with 3 values) with 100,-200,100. "yours" has always 100. You're matrix has 30,30 size, which is normal, because it is only the first part of cycle.
I forgot to mention @Oli in last comment =)
|

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.