make the execution faster

8 views (last 30 days)
prabal datta
prabal datta on 19 Oct 2024
Edited: Bruno Luong on 21 Oct 2024
I wrote the follwoing code . But it is taking too much time to execute, particularly the for loop. Is it possible to reduce the execution time by vectorizing the for loop?
clear all;
close all;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%
x0=1;y0=1;z0=1;
h=0.001;
tend=100000;
tdis= 90000;
t=0:h:tend;
p=round(tdis/h);
x=zeros(3, length(t));
x(:, 1)=[x0; y0; z0];
for i=1:length(t)-1
k1=func(x(:, i));
k2=func(x(:,i)+0.5*h*k1);
k3=func(x(:,i)+0.5*h*k2);
k4=func(x(:,i)+h*k3);
x(:,i+1)=x(:,i)+(1/6)*(k1+2*(k2+k3)+k4)*h;
end
sol_x=x(1, p:end);
sol_z=x(3, p:end);
plot(sol_x, sol_z, 'b', 'Linewidth',0.5);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xdot=func(x)
sigma=10;
b=8/3;
r=28;
xdot=[sigma*(x(2,:)-x(1,:));
r*x(1,:)-x(2,:)-x(1,:)*x(3,:);
x(1,:)*x(2,:)-b*x(3,:)
];
end
  3 Comments
Bruno Luong
Bruno Luong on 20 Oct 2024
Edited: Bruno Luong on 20 Oct 2024
Obviously this is Lorenz ode system using runge kutta 4th order numerial scheme.
埃博拉酱
埃博拉酱 on 21 Oct 2024
There is no way to avoid the computational resource consumption of high-precision iterations for general ODE problems. All I can think of is to use MATLAB built-in functions wherever possible, or to call C++ libraries using MEX file functions, and even then you can't expect very significant performance gains.

Sign in to comment.

Answers (4)

Sam Chak
Sam Chak on 20 Oct 2024
You might as well use the built-in ODE solver, such as ode45, to solve the chaotic Lorenz system for comparison of results purposes. @John D'Errico, now I recall assigning a value to a variable named beta, even though there is a special Beta function in MATLAB.
%% use built-in ODE solver
tspan = [0 150];
x0 = [1; 1; 1];
[t, x] = ode45(@Lorenz, tspan, x0);
%% plot result
plot3(x(:,1), x(:,2), x(:,3)), grid on
xlabel('x_1'), ylabel('x_2'), zlabel('x_3')
view(45, 30)
%% The Chaotic Lorenz System
function xdot = Lorenz(t, x)
sigma = 10;
beta = 8/3;
rho = 28;
xdot = [sigma*(x(2,:) - x(1,:));
rho*x(1,:) - x(2,:) - x(1,:)*x(3,:);
x(1,:)*x(2,:) - beta*x(3,:)];
end
  3 Comments
Bruno Luong
Bruno Luong on 20 Oct 2024
Edited: Bruno Luong on 20 Oct 2024
Observe how the two numerical sollutions diverge after t >= 13
close all
x0=1;y0=1;z0=1;
h=0.0001;
tend=100; % <--- is 150 good enough, instead of 100,000?
tmanual=0:h:tend;
x=zeros(3, length(tmanual));
x(:, 1)=[x0; y0; z0];
for i=1:length(tmanual)-1
k1=func(x(:, i));
k2=func(x(:,i)+0.5*h*k1);
k3=func(x(:,i)+0.5*h*k2);
k4=func(x(:,i)+h*k3);
x(:,i+1)=x(:,i)+(1/6)*(k1+2*(k2+k3)+k4)*h;
end
xmanual = x';
%% use built-in ODE solver
tspan = [0 tend];
xyz0 = [x0; y0; z0];
[t, xyzmatlab] = ode45(@Lorenz, tspan, xyz0);
% interpolation to make the same time basis
xmanual = interp1(tmanual, xmanual, t, 'spline', 'extrap');
% distance of two numerical solutions
dxyz = vecnorm(xmanual-xyzmatlab,2,2);
%% animate result
figure
set(gcf, 'Position', [100 100 1000 430])
subplot(1,2,1)
plot3(xyzmatlab(:,1), xyzmatlab(:,2), xyzmatlab(:,3), 'Color', 0.7+[0 0 0]);
hold on
view(13,10)
axis([ -20 20 -25 30 0 50])
xlabel('x'), ylabel('y'), zlabel('z')
an1 = animatedline(xyzmatlab(1,1), xyzmatlab(1,2), xyzmatlab(1,3), 'MaximumNumPoints', 100, 'Color', 'b','Linewidth', 2);
an2 = animatedline(xmanual(1,1), xmanual(1,2), xmanual(1,3), 'MaximumNumPoints', 100, 'Color', 'r', 'Linewidth', 2);
drawnow
subplot(1,2,2)
axis([0 tend 0 60])
xlabel('t')
ylabel('dxyz')
an3 = animatedline(t(1),dxyz(1), 'Color', 'k');
for i=2:length(t)
addpoints(an1, xyzmatlab(i,1), xyzmatlab(i,2), xyzmatlab(i,3))
addpoints(an2, xmanual(i,1), xmanual(i,2), xmanual(i,3))
addpoints(an3, t(i),dxyz(i))
pause(0.02)
drawnow
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xdot=func(x)
% The Chaotic Lorenz System
sigma=10;
b=8/3;
r=28;
xdot=[sigma*(x(2,:)-x(1,:));
r*x(1,:)-x(2,:)-x(1,:)*x(3,:);
x(1,:)*x(2,:)-b*x(3,:)
];
end
function xdot = Lorenz(t, x)
xdot=func(x);
end
prabal datta
prabal datta on 21 Oct 2024
Edited: prabal datta on 21 Oct 2024
Thank you for the suggession. Unfortunately, I cannot decrease tend (tend=100000), as I need to use the time series for further calculation. tend=100 or 150 is not sufficient for my work. Any further suggessions to make my code faster without affecting tend, h and the methodology (RK4). Checked that most of the time is taking for the execution of the 'for' loop.

Sign in to comment.


Rik
Rik on 19 Oct 2024
There isn't something specific I notice in your loop code (other than your abundance of colon-indexing).
So let's take a look at the rest of your code:
h=0.001;
tend=100000;
t=0:h:tend;
iterations = numel(t)
iterations = 100000001
That seems a lot of iterations. If you want to run your code in a minute, how much time would each iteration be allowed to last?
milisec_per_it = 60*1000/iterations
milisec_per_it = 6.0000e-04
Oof. Less than a 1000th of a milisecond. That seems very short.
Let's park this for now and look at the end:
tdis= 90000;
p=round(tdis/h);
x=rand(3, length(t)); % fake some data
sol_x=x(1, p:end);
sol_z=x(3, p:end);
plotted_points = numel(sol_x)
plotted_points = 10000002
I can't read this easily, so let's convert to exponential:
fprintf('%.1e\n',plotted_points)
1.0e+07
You're plotting 10 million points. A 4k screen has 8.3 million pixels. So unless you have you have monstrous screens with extreme resolutions, you are not going to be able to distinguis any of these points, let along the lines drawn between them.
In conclusion:
You're calculating an enourmous number of iterations. This just takes time. You simply can't finish 10 million calculations in 5 miliseconds. You are also plotting a ridiculous number of points.
Consider using linspace and playing with your step size to reduce the number of iterations and the number of plotted points. You might also want to consider plotting the points as points, instead of the default lines.
  2 Comments
prabal datta
prabal datta on 19 Oct 2024
We can skip the plotting. Any other suggessions?
Rik
Rik on 20 Oct 2024
The fundamental problem is that you want to do a huge amount of calculations. There are some small things you can do (as Walter is helping you do), but the fundamental thing is that this will keep taking a long time to do. It is unclear to me that you actually realize this.

Sign in to comment.


Walter Roberson
Walter Roberson on 19 Oct 2024
Is it possible to reduce the execution time by vectorizing the for loop?
for i=1:length(t)-1
k1=func(x(:, i));
k2=func(x(:,i)+0.5*h*k1);
k3=func(x(:,i)+0.5*h*k2);
k4=func(x(:,i)+h*k3);
x(:,i+1)=x(:,i)+(1/6)*(k1+2*(k2+k3)+k4)*h;
end
The input value, x(:,i) of one step is dependent on the value calculated on the previous iteration. That is not something that can be vectorized.
The most you could do would be to "unroll" the loop -- have your code calculate two (or more) iterations at one time, to reduce the overhead of the for loop.
  4 Comments
prabal datta
prabal datta on 20 Oct 2024
Thank you. But, unfortunately, the "unroll" code is taking almost same time. Any other suggession?
Walter Roberson
Walter Roberson on 20 Oct 2024
Edited: Walter Roberson on 20 Oct 2024
Your possibilities after that involve vectorizing func() and calculating several coefficients at the same time.
Actually it looks like func() is already vectorized, so
k14 = func([x(:, i), x(:,i)+0.5*h*k1, x(:,i)+0.5*h*k2, x(:,i)+h*k3);
x(:,i+1) = x(:,i)+(1/6)*(k14(:,1)+2*(k14(:,2)+k14(:,3))+k14(:,4))*h;

Sign in to comment.


Bruno Luong
Bruno Luong on 21 Oct 2024
Edited: Bruno Luong on 21 Oct 2024
Using ode45 is a better choice because it adapts the time step intelligently.
On my laptop it lasts less the 4 seconds for tend=100000 to solve the ode.
% Elapsed time is 3.858248 seconds.
close all
x0=1;y0=1;z0=1;
tend=100000; % <--- is 150 good enough, instead of 100,000?
sigma=10;
b=8/3;
r=28;
%% use built-in ODE solver
tspan = [0 tend];
xyz0 = [x0; y0; z0];
tic
[t, xyzmatlab] = ode45(@(t,x) func(x, sigma, b, r), tspan, xyz0);
toc
Elapsed time is 49.176212 seconds.
% interpolation to regular time span
h=0.001;
ti = 0:h:tend;
x = interp1(t, xyzmatlab, ti, 'spline', 'extrap');
%% plot result
figure
plot3(x(:,1), x(:,2), x(:,3), 'Color', 0.7+[0 0 0]);
view(13,10)
axis([ -20 20 -25 30 0 50])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xdot=func(x, sigma, b, r)
xdot=[sigma*(x(2)-x(1));
r*x(1)-x(2)-x(1)*x(3);
x(1)*x(2)-b*x(3)
];
end
  1 Comment
Bruno Luong
Bruno Luong on 21 Oct 2024
Edited: Bruno Luong on 21 Oct 2024
The problem is that you (or your advisor) wantt to compute an unneccesary very dense data on a long period of time. If you want to do something quick start with a smart specification with the problem.

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!