Matlab使用ODE45进行积分的向量化处理，包括if和for循环。

huangapple
go评论81阅读模式
英文:

Matlab vectorization with if and for loops using ODE45 to integrate

问题

I am interested in optimizing the speed of my code and when using "Run and Time" this is the function in my code that greatly impacts speed, but I have a hard time conceptualizing how to vectorize this function properly as I usually just do looping, in the attempt I've made I run into an error as it is also used in an integration, my original function is as follows and does not result in an error

  1. function [dotStates] = ODEFunc(t,states,params)
  2. %ODE function
  3. % Loading in and assigning the variables from parameters
  4. K = params(1);                                                             
  5. N = params(2);
  6. nn = params(3);
  7. % magnitude of the coupling based on the number of neighbours
  8. kn = K/nn;
  9. w = params(4:end);
  10. dotStates=states;
  11.     % For each oscillator
  12.     for i=1:N
  13.         % Use the oscillators natural frequency
  14.         dotStates(i) = w(i);
  15.         % For j number of neighbours
  16.         for j=(i-nn):(i+nn)
  17.             % neighbour number is positive and shorter than # of oscilators
  18.             if (j > 0) && (j < length(dotStates))
  19.                 dotStates(i) =  dotStates(i) + (kn * sin( states(j)-states(i) ));
  20.             end
  21.         end
  22.     end
  23. end

I've tried following the mathworks vectorization guide: MathWorks Vectorization Guide

My attempt so far has been to follow some of the inputs of what they use, such as using a mask and have generated the following code

  1. function [dotStates] = ODEFunc(t,states,params)
  2. %ODE function
  3. % Loading in and assigning the variables from parameters
  4. K = params(1);                                                             
  5. N = params(2);
  6. nn = params(3);
  7. % magnitude of the coupling based on the number of neighbours
  8. kn = K/nn;
  9. w = params(4:end);
  10. dotStates=states;
  11. % Use the oscillators natural frequency
  12. dotStates = w';
  13. % Mask of j states
  14. j = (i-nn):(i+nn);
  15. % neighbours cannot exceed boundaries
  16. j = j(j>0 & j <= length(dotStates));
  17. jstate = states(j);
  18. jstate(numel(states)) = 0;
  19. dotStates =  dotStates + (kn * sin( jstate'-states ));
  20. end

I ended up with a vector that is shorter than what is being written to, and my solution has been to just add a bunch of zeros to the "jstate" variable to make up for the difference, but that does not feel like proper vectorization. When I run the code, I get the following error which is tied to an integration step:

Warning: Colon operands must be real scalars.
In RK_ODE_2411>ODEFunc (line 99)
In RK_ODE_2411>@(t,states)ODEFunc(t,states,params)
In ode45 (line 324)
In RK_ODE_2411 (line 58)

The function is in turn used in the following segment for the integration using ODE45

  1. %% Integration via ODE45
  2. for K = 0:.1:Klen
  3.     params(1) = K;
  4.     K_count = K_count+1;
  5.     nn_count = 0;
  6.     for nn = nnlen:nnlen
  7.         params(3) = nn;
  8.         % index counter
  9.         nn_count = nn_count+1;
  10.         % 6th order runge kutta
  12.         sol(K_count,nn_count) = ode45(@(t,states) ODEFunc(t,states,params),tSpan,init,options);
  13.     end
  14. end

Where line 58 is

  1. sol(K_count,nn_count) = ode45(@(t,states) ODEFunc(t,states,params),tSpan,init,options);

EDIT: Line 99 in ODEFunc is

  1. j = (i-nn):(i+nn);
英文:

I am interested in optimizing the speed of my code and when using "Run and Time" this is the function in my code that greatly impacts speed, but I have a hard time conceptualizing how to vectorize this function properly as I usually just do looping, in the attempt I've made I run into an error as it is also used in an integration, my original function is as follows and does not result in an error

  1. function [dotStates] = ODEFunc(t,states,params)
  2. %ODE function
  3. % Loading in and assigning the variables from parameters
  4. K = params(1);                                                             
  5. N = params(2);
  6. nn = params(3);
  7. % magnitude of the coupling based on the number of neighbours
  8. kn = K/nn;
  9. w = params(4:end);
  10. dotStates=states;
  11.     % For each oscillator
  12.     for i=1:N
  13.         % Use the oscillators natural frequency
  14.         dotStates(i) = w(i);
  15.         % For j number of neighbours
  16.         for j=(i-nn):(i+nn)
  17.             % neighbour number is positive and shorter than # of oscilators
  18.             if (j &gt; 0) &amp;&amp; (j &lt; length(dotStates))
  19.                 dotStates(i) =  dotStates(i) + (kn * sin( states(j)-states(i) ));
  20.             end
  21.         end
  22.     end
  23. end

I've tried following the mathworks vectorization guide: https://se.mathworks.com/help/matlab/matlab_prog/vectorization.html

My attempt so far has been to follow some of the inputs of what they use, such as using a mask and have generated following code

  1. function [dotStates] = ODEFunc(t,states,params)
  2. %ODE function
  3. % Loading in and assigning the variables from parameters
  4. K = params(1);                                                             
  5. N = params(2);
  6. nn = params(3);
  7. % magnitude of the coupling based on the number of neighbours
  8. kn = K/nn;
  9. w = params(4:end);
  10. dotStates=states;
  11. % Use the oscillators natural frequency
  12. dotStates = w&#39;;
  13. % Mask of j states
  14. j = (i-nn):(i+nn);
  15. % neighbours cannot exceed boundaries
  16. j = j(j&gt;0 &amp; j &lt;=length(dotStates));
  17. jstate = states(j);
  18. jstate(numel(states)) = 0;
  19. dotStates =  dotStates + (kn * sin( jstate&#39;-states ));
  20. end

I ended up with a vector that is shorter than what is being written to and my solution has been to just add a bunch of zeros to the "jstate" variable to make up for the difference but that does not feel like a proper vectorization and when I run the code I get the following error which is tied to and integration step

>Warning: Colon operands must be real scalars.

> In RK_ODE_2411>ODEFunc (line 99)
In RK_ODE_2411>@(t,states)ODEFunc(t,states,params)
In ode45 (line 324)
In RK_ODE_2411 (line 58)

the function is in turn used in the following segment for the integration using ODE45

  1. %% Integration via ODE45
  2. for K = 0:.1:Klen
  3.     params(1) = K;
  4.     K_count = K_count+1;
  5.     nn_count = 0;
  6.     for nn = nnlen:nnlen
  7.         params(3) = nn;
  8.         % index counter
  9.         nn_count = nn_count+1;
  10.         % 6th order runge kutta
  11.         
  12.         sol(K_count,nn_count) = ode45(@(t,states) ODEFunc(t,states,params),tSpan,init,options);
  13.     end
  14. end

where line 58 is

  1. sol(K_count,nn_count) = ode45(@(t,states) ODEFunc(t,states,params),tSpan,init,options);

EDIT: line 99 in ODEFunc is

  1. j = (i-nn):(i+nn);

答案1

得分: 1

尝试这段代码片段

  1. % 对于每个振荡器
  2. for i = 1:N
  3.     % 对于邻居数为j
  4.     j = (i - nn):(i + nn);
  6.     % 邻居数是正数且比振荡器的数量短
  7.     lg = (j > 0) & (j < length(dotStates));
  8.     dotStates(i) = w(i) + sum(kn * sin(states(lg) - states(i)));
  9. end

最重要的是，确保 dotStates 不会大于 stats，因为这会迫使Matlab重新排列内存，这会严重减慢代码的运行速度。

英文:

Try this snippet

  1. % For each oscillator
  2. for i=1:N
  3.     % For j number of neighbours
  4.     j=(i-nn):(i+nn);
  5.     
  6.     % neighbour number is positive and shorter than # of oscilators
  7.     lg = (j &gt; 0) &amp; (j &lt; length(dotStates));
  8.     dotStates(i) =  w(i) + sum(kn * sin( states(lg)-states(i) ));
  9. end

the most important is though that dotStates won't be larger than stats, since this would force matlab to rearrange its memory, which slows down the code enormously.

huangapple
  • 本文由 发表于 2020年1月6日 18:49:04
  • 转载请务必保留本文链接：https://go.coder-hub.com/59610737.html
  • if-statement
  • matlab
  • nested-loops
  • ode45
  • vectorization
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定