# Improving MATLAB programme to make it Run faster

2 views (last 30 days)
aakash dewangan on 12 Nov 2021
Commented: aakash dewangan on 15 Nov 2021
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 on 12 Nov 2021

Jan on 12 Nov 2021
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 on 15 Nov 2021
Thank you very much for your time and help :)

### Categories

Find more on Argument Definitions in Help Center and File Exchange

R2021a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!