nested loops including ode45: make code faster?
in [Matlab]
|
From: Alexander Erlich on 19 Jan 2010 05:55 Hi, I have written a program which uses nested loops in which an ode45 command is executed. Then, a dot of a certain color is plotted. I have attatched the source codes and also uploaded them as m-files (along with a picture of the time consuming final plot) here: http://www.airlich.de/filebase/nonlinear_dynamics/Duffing_attractors/FEX_question/ My problem is, this code is very slow, and I would like to make it faster. The idea is this: I have two nested for loops, x=-4:0.1:4 and y=-2:0.1:2. Within every y loop, I call ode45 with [x;y] als the initial conditions. The equation which I am investigating here is the Duffing oscillator equation. After ode45 has computed all the trajectories, I only take the final point of one trajectory (namely, x1. So I compute x1(end)) because I am only interested in the long term behaviour of the equation (possibly, there are more efficient ways than computing the whole solution just to find the behaviour after, say, 100sec?). Finally, I check whether the computed x1(end) value is close to the attractor x1--->1 or x1--->-1. In the first case, I plot a blue dot at (x,y), in the latter case, a yellow dot. My feeling tells me that in this program, not the plotting or the loops are the delimiting factor, but ode45 (which is already a very fast algorithm) is. But I am not quite sure. I'd be very grateful if you could tell me how I can make the code faster, because I believe that for my Bachelor thesis I am going to write many more similar programs in order to study dynamical systems, so it would be good to know how to make the code reasonably fast. Alexander ----------------------------------Duffing_attractors------------------------------- function Duffing_attractors(interval, gamma, tol) %DUFFING_attractors Analyzes the attractor basins of the Duffing oscillator. % % Example call: Duffing_attractors([0, 100],0,0.1); % % The corresponding ODEs are in the file DuffingDGL.m. They are taken from % Argyris / Faust / Haase: Die Erforschung des Chaos. Vieweg, 1994, p. 669. % The parameters (C, beta, alpha, gamma, omega_E) in the example call are % chosen to be the same as in Argyris, p. 673. global C beta alpha gamma omega_E; C=0.12; beta=0.5; alpha=0.5; omega_E=0.7; clc; clf; clear All; % clear command, figure, all variables figure(1); h=plot(1,1); hold on; for x=-4:0.1:4 if(ishandle(h)==0) break; end for y=-2:0.1:2 if(ishandle(h)==0) % checks if figure has been closed. If so, ... break; % ... simply leave the loop end [t,Xcurrent]=ode45(@DuffingODE,interval, [x;y]); x1=Xcurrent(:,1); % x2=Xcurrent(:,2); % x2 = dx1/dt will not be required. % if long-term behaviour is attracted to 1 (i.e. [1,0] in the % phase space portrait), draw a blue dot. if(abs(x1(end)-1)<tol) plot(x,y,'o','Color','Blue','linewidth',6) drawnow; %update the plot end % if long-term behaviour is attracted to -1 (i.e. [-1,0] in the % phase space portrait), draw a yellow dot. if(abs(x1(end)+1)<tol) plot(x,y,'o','Color','Yellow','linewidth',6) drawnow; %update the plot end end end ----------------------------------DuffingODE------------------------------- function xdot = DuffingODE(t,x) %DuffingODE contains set of two ODEs for Duffing oscillator. global C beta alpha gamma omega_E; xdot=zeros(2,1); xdot(1)=x(2); xdot(2)=-C*x(2)+beta*x(1)-alpha*x(1)^3+gamma*cos(omega_E*t);
From: Alexander Erlich on 26 Jan 2010 18:50 If you try e.g. >> Duffing_attractors([0, 1],0,0.1); (this means you solve each ode45(...) for t=1 instead of t=100), the program is still remarkably slow. What makes it so slow? Is it ode45, the nested loops, the plotting (hold on, drawnow...), or maybe something else? What would you suggest? "Alexander Erlich" <alexander.erlich(a)gmail.com> wrote in message <hj4328$mde$1(a)fred.mathworks.com>... > Hi, > > I have written a program which uses nested loops in which an ode45 command is executed. Then, a dot of a certain color is plotted. I have attatched the source codes and also uploaded them as m-files (along with a picture of the time consuming final plot) here: > > http://www.airlich.de/filebase/nonlinear_dynamics/Duffing_attractors/FEX_question/ > > My problem is, this code is very slow, and I would like to make it faster. The idea is this: I have two nested for loops, x=-4:0.1:4 and y=-2:0.1:2. Within every y loop, I call ode45 with [x;y] als the initial conditions. The equation which I am investigating here is the Duffing oscillator equation. After ode45 has computed all the trajectories, I only take the final point of one trajectory (namely, x1. So I compute x1(end)) because I am only interested in the long term behaviour of the equation (possibly, there are more efficient ways than computing the whole solution just to find the behaviour after, say, 100sec?). Finally, I check whether the computed x1(end) value is close to the attractor x1--->1 or x1--->-1. In the first case, I plot a blue dot at (x,y), in the latter case, a yellow dot. > > My feeling tells me that in this program, not the plotting or the loops are the delimiting factor, but ode45 (which is already a very fast algorithm) is. But I am not quite sure. I'd be very grateful if you could tell me how I can make the code faster, because I believe that for my Bachelor thesis I am going to write many more similar programs in order to study dynamical systems, so it would be good to know how to make the code reasonably fast. > > Alexander > > ----------------------------------Duffing_attractors------------------------------- > function Duffing_attractors(interval, gamma, tol) > %DUFFING_attractors Analyzes the attractor basins of the Duffing oscillator. > % > % Example call: Duffing_attractors([0, 100],0,0.1); > % > % The corresponding ODEs are in the file DuffingDGL.m. They are taken from > % Argyris / Faust / Haase: Die Erforschung des Chaos. Vieweg, 1994, p. 669. > % The parameters (C, beta, alpha, gamma, omega_E) in the example call are > % chosen to be the same as in Argyris, p. 673. > > global C beta alpha gamma omega_E; > C=0.12; beta=0.5; alpha=0.5; omega_E=0.7; > > clc; clf; clear All; % clear command, figure, all variables > > figure(1); > > h=plot(1,1); hold on; > > for x=-4:0.1:4 > > if(ishandle(h)==0) > break; > end > > for y=-2:0.1:2 > > if(ishandle(h)==0) % checks if figure has been closed. If so, ... > > break; % ... simply leave the loop > end > > [t,Xcurrent]=ode45(@DuffingODE,interval, [x;y]); > x1=Xcurrent(:,1); > % x2=Xcurrent(:,2); % x2 = dx1/dt will not be required. > > % if long-term behaviour is attracted to 1 (i.e. [1,0] in the > % phase space portrait), draw a blue dot. > > if(abs(x1(end)-1)<tol) > plot(x,y,'o','Color','Blue','linewidth',6) > drawnow; %update the plot > end > > % if long-term behaviour is attracted to -1 (i.e. [-1,0] in the > % phase space portrait), draw a yellow dot. > > if(abs(x1(end)+1)<tol) > plot(x,y,'o','Color','Yellow','linewidth',6) > drawnow; %update the plot > end > end > end > > ----------------------------------DuffingODE------------------------------- > function xdot = DuffingODE(t,x) > %DuffingODE contains set of two ODEs for Duffing oscillator. > > global C beta alpha gamma omega_E; > > xdot=zeros(2,1); > > xdot(1)=x(2); > xdot(2)=-C*x(2)+beta*x(1)-alpha*x(1)^3+gamma*cos(omega_E*t);
From: Jan Simon on 26 Jan 2010 19:06 Dear Alexander! > If you try e.g. > > >> Duffing_attractors([0, 1],0,0.1); > > (this means you solve each ode45(...) for t=1 instead of t=100), the program is still remarkably slow. What makes it so slow? Do you know the profiler? profile on Duffing_attractors([0, 1],0,0.1); profile report And now you can look where the time is vanishing. My opinion: ODE45 is not fast. Use an integrator written in C. Good luck, Jan
From: Jan Simon on 26 Jan 2010 19:38 Dear Alexander - again! > ----------------------------------DuffingODE------------------------------- > function xdot = DuffingODE(t,x) > global C beta alpha gamma omega_E; > > xdot=zeros(2,1); > xdot(1)=x(2); > xdot(2)=-C*x(2)+beta*x(1)-alpha*x(1)^3+gamma*cos(omega_E*t); Save some time here! 1. Do not use global variables. Try if it is faster use further parameters to your function: function xdot = DuffingODE(t,x, C, alpha, beta, gamma, omega_E) And append the parameters to the ODE45 call. 2. Omit the dependency to omega_E. Integrating from 0 to 100 over cos(omega_E*t) means integrating from 0 to 100*omega_E over cos(t). Just the result differs by the factor omega_E. --- I do not think that this single multiplication in the ODEFCN results in a noticiably acceleration. 3. Look on the stepsizes of the integration: Does the estimation of the initial step size work well? I do not expect Matlab's ODE45 to cause troubles with a bad initial guess. But this is often a problem with handmade integrators. On demand define a better stepsize with ODESET. 4. Vectorize the ODEFCN: See "help odeset" -> Vectorized. This means a remarkable speedup in general. Good luck, Jan
From: Alexander Erlich on 4 Feb 2010 09:27
Thank you Jan! I believe I know how to fix 1. and 2. But concerning esp. 4: if it is not too comprehensive/does not involve much research, can you give me some more instructions what to do here? I took a look at the Vectorized section of the odeset help page, but actually didn't quite understand what to do and how to set that vectorized property. Can you give me some more instructions here? Another thing: You recommended fast C integrator. However, I am much more fluent in FORTRAN as I recently did a project in this language. Do you happen to know a FORTRAN integrator which would be very fast? Thank you in advance! Alexander |