make the execution faster

21 次查看（过去 30 天）
prabal datta
prabal datta 2024-10-19
编辑： Bruno Luong 2024-10-21，8:26
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 个评论
显示 1更早的评论
Bruno Luong
Bruno Luong 2024-10-20
编辑：Bruno Luong 2024-10-20
Obviously this is Lorenz ode system using runge kutta 4th order numerial scheme.
埃博拉酱
埃博拉酱 2024-10-21，2:46
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（4 个）

Sam Chak
Sam Chak 2024-10-20
Ran in:
Hi @prabal datta
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 个评论
显示 1更早的评论
Bruno Luong
Bruno Luong 2024-10-20
编辑：Bruno Luong 2024-10-20
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 2024-10-21，2:31
编辑：prabal datta 2024-10-21，2:36
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.

请先登录，再进行评论。

Rik
Rik 2024-10-19
Ran in:
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 个评论
prabal datta
prabal datta 2024-10-19
We can skip the plotting. Any other suggessions?
Rik
Rik 2024-10-20
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.

请先登录，再进行评论。

Walter Roberson
Walter Roberson 2024-10-19
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 个评论
显示 2更早的评论
prabal datta
prabal datta 2024-10-20
Thank you. But, unfortunately, the "unroll" code is taking almost same time. Any other suggession?
Walter Roberson
Walter Roberson 2024-10-20
编辑：Walter Roberson 2024-10-20
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;

请先登录，再进行评论。

Bruno Luong
Bruno Luong 2024-10-21，7:32
编辑：Bruno Luong 2024-10-21，8:19
Ran in:
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 个评论
Bruno Luong
Bruno Luong 2024-10-21，8:25
编辑：Bruno Luong 2024-10-21，8:26
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by