make the execution faster
8 views (last 30 days)
Show older comments
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
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.
Answers (4)
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
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
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)
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
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)
I can't read this easily, so let's convert to exponential:
fprintf('%.1e\n',plotted_points)
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
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.
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
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;
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
% 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
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.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!