make the execution faster

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  

[Steven Lord]

Why are you taking such tiny steps over such a long time period (relatively speaking)? What about the problem you're trying to solve requires it?

h=0.001;
tend=100000;

If this isn't for a homework assignment, why not use the functions included in MATLAB to solve this type of problem?

[Bruno Luong]

Obviously this is Lorenz ode system using runge kutta 4th order numerial scheme.

[埃博拉酱]

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.

Answer 1

[Sam Chak]

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  

[Sam Chak]

Hi @prabal datta

Is tend = 150 good enough, instead of 100,000?

x0=1;y0=1;z0=1;
h=0.001;
tend=150;           % <--- is 150 good enough, instead of 100,000?
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, 1:end);  % <--- start from 1, instead of p
sol_z=x(3, 1:end);  % <--- start from 1, instead of p
plot(sol_x, sol_z, 'b', 'Linewidth',0.5), grid on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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  

[Bruno Luong]

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]

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.

Answer 2

[Rik]

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.

[prabal datta]

We can skip the plotting. Any other suggessions?

[Rik]

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.

Answer 3

[Walter Roberson]

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.

[prabal datta]

Thank you for the suggession. How can I "unroll" the loop or any other modification to speed-up the code?

[Walter Roberson]

for i=1:2: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;
    k1=func(x(:, i+1));
    k2=func(x(:,i+1)+0.5*h*k1);
    k3=func(x(:,i+1)+0.5*h*k2);
    k4=func(x(:,i+1)+h*k3);
    x(:,i+2)=x(:,i+1)+(1/6)*(k1+2*(k2+k3)+k4)*h;
end

It is unlikely that this will save significant time. Also, it needs a trailing section to potentially do the last iteration depending on whether length(t) is even or odd

[prabal datta]

Thank you. But, unfortunately, the "unroll" code is taking almost same time. Any other suggession?

[Walter Roberson]

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;

Answer 4

[Bruno Luong]

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  

[Bruno Luong]

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.