accelerating code in matlab - for loop

Asked by Sara on 5 Jul 2012
Latest activity Commented on by Matt Kindig on 12 Jul 2012

Hi Guys,

    int = zeros(1,length(x));
    if k <= 0
        int = x ;
        return
    end
    for n = 3:length(x)
          y1 = x(2:n) ;
          t1 = ((n-2):-1:0)*dt ;
         y2 = x(1:n-1) ;
          t2 = t1 +dt; 
          int(n) = sum (t1.^(k-1)/factorial(k-1).*y1 + t2.^(k-1)/factorial(k-1).*y2)*dt/2 ;
    end

is there any suggestion to write this part of code more optimal? It is used for computing Integration. The most time consuming line is int(n)= ... Whether I can change and remove the for loop here, and implement the for functionality in another way, e.g. vectors? Any suggestions would be appreciated ...

2 Comments

Ryan on 5 Jul 2012
for n = 3:length(x ) 
    y1 = x(2:n) ; 
    t1 = ((n-2):-1:0)*dt ; 
    y2 = x(1:n-1) ; 
    t2 = t1 +dt; 
    int(n) = sum (t1.^(k-1)/factorial(k-1).*y1 + t2.^(k-1)/factorial(k-...
                 1).*y2)*dt/2 ;
end
Jan Simon on 11 Jul 2012

@Sara: You can delete your question using the "Delete" button left beside the question. If such a button does not appear although you are logged in, please contact an editor or files@mathworks.com and ask for deleteing the thread. Duplicate posts waste the time of the ones, who want to help you.

Sara

Products

No products are associated with this question.

2 Answers

Answer by Jan Simon on 11 Jul 2012
Accepted answer

At first I've clean up your code a bit:

if k <= 0
   int = x ;
   return
end
int = zeros(1, length(x));
fk1 = factorial(k - 1);
dt2 = dt / 2;
for n = 3:length(x)
  y1 = x(2:n);
  t1 = ((n-2):-1:0) * dt;
  y2 = x(1:n-1);
  t2 = t1 + dt; 
  int(n) = sum(t1 .^ (k-1) / fk1 .* y1 + t2 .^ (k-1) / fk1 .* y2) * dt2;
end

Now some improvements:

lenx = length(x);
fk1 = factorial(k - 1);
t1Vec = (((lenx - 2):-1:0) * dt) .^ (k-1) * dt / fk1;
int = zeros(1, lenx);
dt2 = dt / 2;
for n = 3:lenx
  y1 = x(2:n);
  t1 = t1vec(lenx-a-n:lenx-b);  % **see comment**
  y2 = x(1:n-1);
  t2 = t1 + dt; 
  int(n) = sum(t1 .* y1 + t2 .^ (k-1) / fk1 .* y2) * dt2;
end

I do not have the time to find the right constantd for a and b. Without access to Matlab I cannot simply try it, but you can. The idea is to avoid the repeated expensive calculation of t1 .^(k-1), when all elements except for one have been processed already. The same works for t2. I leave it up to you to elaborate this.

3 Comments

Sara on 11 Jul 2012

to find the right constantd for a and b. :: could you let me know which role these a and b have ? and how can I find the right one of them?

Sara on 11 Jul 2012
t1Vec = (((lenx - 2):-1:0) * dt) .^ (k-1) * dt / fk1;

the second dt is required where .. (k-1)times dt...

Jan Simon on 11 Jul 2012

I assume this is correct, Sara. "a" and "b" are used only to replicate the limits of t1 = ((n-2):-1:0) * dt.

Jan Simon
Answer by Matt Kindig on 5 Jul 2012

Could you please format your code properly? It is difficult to tell which lines are commented, etc. Please edit your question to include proper line breaks, and then apply the 'code' tag (see the {} Code icon in the editor window).

Also, as a first observation, you should pre-allocate your variable 'int'--that is the reason that it is taking so long. Prior to the loop, add this line:

int = NaN(length(x), 1); 

I think that this change alone will substantially improve your performance.

2 Comments

Jan Simon on 11 Jul 2012

int is pre-allocated by zeros already.

Matt Kindig on 12 Jul 2012

Ah yes, I think I missed that before.

Matt Kindig

Contact us