Improving MATLAB programme to make it Run faster

Hello all,

I have written a MATLAB programme to plot my simulation results. Programme is as follows:

tic

clc; clear all; close all

% ==============================================================

axis tight manual

vid = VideoWriter('3D Stability My EOM.avi')

open(vid);

% ========================== ====================================

syms p1(t) p2(t) p3(t) alpha_by_L beta_by_L gamma_by_L epsilon_by_L

Dp1 = diff(p1); D2p1 = diff(p1,2); Dp2 = diff(p2); D2p2 = diff(p2,2); Dp3 = diff(p3); D2p3 = diff(p3,2);

for gamma_by_L = 0.01:0.1:5

for beta_by_L = 0.1:0.1:5

for alpha_by_L = 0:0.1:5

% matrix terms

AA = 1/2 + (alpha_by_L)*(sin(pi*beta_by_L*t))^2;

BB = (alpha_by_L)*sin(pi*beta_by_L*t)*sin(2*pi*beta_by_L*t);

CC = (alpha_by_L)*sin(pi*beta_by_L*t)*sin(3*pi*beta_by_L*t);

DD = 1/2 + (alpha_by_L)*(sin(2*pi*beta_by_L*t))^2;

EE = (alpha_by_L)*sin(3*pi*beta_by_L*t)*sin(2*pi*beta_by_L*t);

FF = 1/2 + (alpha_by_L)*(sin(3*pi*beta_by_L*t))^2;

% matrix terms

GG = (pi^2)/2 + (gamma_by_L)*(sin(pi*beta_by_L*t))^2;

HH = (gamma_by_L)*sin(pi*beta_by_L*t)*sin(2*pi*beta_by_L*t);

II = (gamma_by_L)*sin(3*pi*beta_by_L*t)*sin(pi*beta_by_L*t);

JJ = (4*pi^2)/2 + (gamma_by_L)*(sin(2*pi*beta_by_L*t))^2;

KK = (gamma_by_L)*sin(2*pi*beta_by_L*t)*sin(3*pi*beta_by_L*t);

LL = (9*pi^2)/2 + (gamma_by_L)*(sin(3*pi*beta_by_L*t))^2;

% RHS

MM = (epsilon_by_L)*sin(pi*beta_by_L*t);

NN = (epsilon_by_L)*sin(2*pi*beta_by_L*t);

OO = (epsilon_by_L)*sin(3*pi*beta_by_L*t);

Eq1 = AA*diff(p1,t,2) + BB*diff(p2,t,2) + CC*diff(p3,t,2) + GG*p1 + HH*p2 + II*p3 == 0; % Equation 1

Eq2 = BB*diff(p1,t,2) + DD*diff(p2,t,2) + EE*diff(p3,t,2) + HH*p1 + JJ*p2 + KK*p3 == 0; % Equation 2

Eq3 = CC*diff(p1,t,2) + EE*diff(p2,t,2) + FF*diff(p3,t,2) + II*p1 + KK*p2 + LL*p3 == 0; % Equation 3

%

[V,S] = odeToVectorField(Eq1, Eq2, Eq3);

ftotal = matlabFunction(V, 'Vars',{'t','Y'}); % Using readymade MATLAB function to solve using ODE 45

interval = [0 2/beta_by_L]; % Time Interval to run the program

%%

for i = 1:6;

if i == 1

y0 = [1 ;0 ;0 ;0 ;0 ;0 ] ; % First IC

elseif i == 2

y0 = [0 ;1 ;0 ;0 ;0 ;0 ] ; % Second IC

elseif i == 3

y0 = [0 ;0 ;1 ;0 ;0 ;0 ] ;

elseif i == 4

y0 = [0 ;0 ;0 ;1 ;0 ;0 ] ;

elseif i == 5

y0 = [0 ;0 ;0 ;0 ;1 ;0 ] ;

elseif i == 6

y0 = [0 ;0 ;0 ;0 ;0 ;1 ] ;

end

xSol = ode45( @(t,Y)ftotal(t,Y),interval,y0); % Using ODE 45 to solve state space form of ODE

tValues = linspace(interval(1),interval(2),1000); % dividing time interval

x1Values = deval(xSol,tValues,1); % number 1 denotes first solution likewise you can mention 2 ,3 & 4 for the next three solutions

x2Values = deval(xSol,tValues,2); % number 2 denotes second solution likewise you can mention 2 ,3 & 4 for the next three solutions

x3Values = deval(xSol,tValues,3);

x4Values = deval(xSol,tValues,4);

x5Values = deval(xSol,tValues,5);

x6Values = deval(xSol,tValues,6);

% plot(tValues,x1Values)

if i == 1

col_11 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];

col_12 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];

elseif i == 2

col_21 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];

col_22 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];

elseif i == 3

col_31 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];

col_32 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];

elseif i == 4

col_41 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];

col_42 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];

elseif i == 5

col_51 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];

col_52 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];

elseif i == 6

col_61 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];

col_62 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];

end

end

A1 = [col_11 col_21 col_31 col_41 col_51 col_61];

A2 = [col_12 col_22 col_32 col_42 col_52 col_62] ;

A = (inv(A1))*A2

eigval = eig(A);

abseigval = abs(eigval)

iteval = gamma_by_L

%

if (abseigval(1)<=1.010 & abseigval(2)<=1.010 & abseigval(3)<=1.0010 & abseigval(4)<=1.0010 & abseigval(5)<=1.0010 & abseigval(6)<=1.0010)

figure(1)

plot(beta_by_L,alpha_by_L,'k.', 'MarkerSize', 10)

xlabel('beta_by_L')

ylabel('alpha_by_L')

zlabel('gamma_by_L')

title('Stable points')

hold on

else

figure(1)

plot(beta_by_L,alpha_by_L,'r.', 'MarkerSize', 10)

xlabel('beta_by_L')

ylabel('alpha_by_L')

zlabel('gamma_by_L')

title('Untable points')

hold on

end

end

end

hold off

xlim([0 5])

ylim([0 5])

zlim([0 5])

frame = getframe(gcf);

writeVideo(vid,frame)

end

close(vid);

toc

here I am trying to plot different combinations of "alpha_by_L" and "beta_by_L" for different values of "gamma_by_L". It took 3 days and still it is not finished yet. Can somebody please help me to improve this code, so that it will run faster on MATLAB.

Really need help!

Thanks :)

[KSSV]

Read about profile. This will help you to know which part will take time.

https://in.mathworks.com/help/matlab/ref/profile.html

Answer 1

[Jan]

In Matlab online one iteration needs about 3 seconds.

for gamma_by_L = 0.01:0.1:5
  for beta_by_L = 0.1:0.1:5
     for alpha_by_L = 0:0.1:5

These are 127'500 iterations. This means about 4.5 days of processing.

As KSSV suggests, use the profiler to find the bottleneck. There is no obvious point, but some small modifications are useful:

• Omit the clear all but move the code into a function to keep the worksapce clean.
• Replace (inv(A1))*A2 by A1 \ A2.
• Suppress the output to the command window by appending semicolons.
• And just nicer without any speedup:

   % Replace:
   if i == 1 
       y0 = [1 ;0 ;0 ;0 ;0 ;0 ] ; % First IC
    
   elseif i == 2
       y0 = [0 ;1 ;0 ;0 ;0 ;0 ] ; % Second IC
 
   elseif i == 3
       y0 = [0 ;0 ;1 ;0 ;0 ;0 ] ;
   
   elseif i == 4
       y0 = [0 ;0 ;0 ;1 ;0 ;0 ] ;
   
   elseif i == 5
       y0 = [0 ;0 ;0 ;0 ;1 ;0 ] ; 
   
   elseif i == 6
       y0 = [0 ;0 ;0 ;0 ;0 ;1 ] ;     
   end
   
   % by:
   y0    = zero(6, 1);
   y0(i) = 1;   

These points will not cause a relevant speedop. The problem is huge.

[aakash dewangan]

Thank you very much for your time and help :)