Search Discussions

Athanasios Paraskevopoulos
Last activity 2024-8-5

Dynamics of predator-prey model

This project discusses predator-prey system, particularly the Lotka-Volterra equations,which model the interaction between two sprecies: prey and predators. Let's solve the Lotka-Volterra equations numerically and visualize the results.% Define parameters

% Define parameters

alpha = 1.0; % Prey birth rate

beta = 0.1; % Predator success rate

gamma = 1.5; % Predator death rate

delta = 0.075; % Predator reproduction rate

% Define the symbolic variables

syms R W

% Define the equations

eq1 = alpha * R - beta * R * W == 0;

eq2 = delta * R * W - gamma * W == 0;

% Solve the equations

equilibriumPoints = solve([eq1, eq2], [R, W]);

% Extract the equilibrium point values

Req = double(equilibriumPoints.R);

Weq = double(equilibriumPoints.W);

% Display the equilibrium points

equilibriumPointsValues = [Req, Weq]

equilibriumPointsValues = 2x2

0 0 20 10

% Solve the differential equations using ode45

lotkaVolterra = @(t,Y)[alpha*Y(1)-beta*Y(1)*Y(2);

delta*Y(1)*Y(2)-gamma*Y(2)];

% Initial conditions

R0 = 40;

W0 = 9;

Y0 = [R0, W0];

tspan = [0, 100];

% Solve the differential equations

[t, Y] = ode45(lotkaVolterra, tspan, Y0);

% Extract the populations

R = Y(:, 1);

W = Y(:, 2);

% Plot the results

figure;

subplot(2,1,1);

plot(t, R, 'r', 'LineWidth', 1.5);

hold on;

plot(t, W, 'b', 'LineWidth', 1.5);

xlabel('Time (months)');

ylabel('Population');

legend('R', 'W');

grid on;

subplot(2,1,2);

plot(R, W, 'k', 'LineWidth', 1.5);

xlabel('R');

ylabel('W');

grid on;

hold on;

plot(Req, Weq, 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r');

legend('Phase Trajectory', 'Equilibrium Point');

Now, we need to handle a modified version of the Lotka-Volterra equations. These modified equations incorporate logistic growth fot the prey population.

These equations are:

% Define parameters

alpha = 1.0;

K = 100; % Carrying Capacity of the prey population

beta = 0.1;

gamma = 1.5;

delta = 0.075;

% Define the symbolic variables

syms R W

% Define the equations

eq1 = alpha*R*(1 - R/K) - beta*R*W == 0;

eq2 = delta*R*W - gamma*W == 0;

% Solve the equations

equilibriumPoints = solve([eq1, eq2], [R, W]);

% Extract the equilibrium point values

Req = double(equilibriumPoints.R);

Weq = double(equilibriumPoints.W);

% Display the equilibrium points

equilibriumPointsValues = [Req, Weq]

equilibriumPointsValues = 3x2

0 0 20 8 100 0

% Solve the differential equations using ode45

modified_lv = @(t,Y)[alpha*Y(1)*(1-Y(1)/K)-beta*Y(1)*Y(2);

delta*Y(1)*Y(2)-gamma*Y(2)];

% Initial conditions

R0 = 40;

W0 = 9;

Y0 = [R0, W0];

tspan = [0, 100];

% Solve the differential equations

[t, Y] = ode45(modified_lv, tspan, Y0);

% Extract the populations

R = Y(:, 1);

W = Y(:, 2);

% Plot the results

figure;

subplot(2,1,1);

plot(t, R, 'r', 'LineWidth', 1.5);

hold on;

plot(t, W, 'b', 'LineWidth', 1.5);

xlabel('Time (months)');

ylabel('Population');

legend('R', 'W');

grid on;

subplot(2,1,2);

plot(R, W, 'k', 'LineWidth', 1.5);

xlabel('R');

ylabel('W');

grid on;

hold on;

plot(Req, Weq, 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r');

legend('Phase Trajectory', 'Equilibrium Point');

Image Analyst
Last activity 2024-8-10

Setting the Microsoft Azure Information Protection label

Does your company or organization require that all your Word Documents and Excel workbooks be labeled with a Microsoft Azure Information Protection label or else they can't be saved? These are the labels that are right below the tool ribbon that apply a category label such as "Public", "Business Use", or "Highly Restricted". If so, you can either

Create and save a "template file" with the desired label and then call copyfile to make a copy of that file and then write your results to the new copy, or
If using Windows you can create and/or open the file using ActiveX and then apply the desired label from your MATLAB program's code.

For #1 you can do

copyfile(templateFileName, newDataFileName);
writematrix(myData, newDataFileName);

If the template has the AIP label applied to it, then the copy will also inherit the same label.

For #2, here is a demo for how to apply the code using ActiveX.

% Test to set the Microsoft Azure Information Protection label on an Excel workbook.
% Reference support article:
% https://www.mathworks.com/matlabcentral/answers/1901140-why-does-azure-information-protection-popup-pause-the-matlab-script-when-i-use-actxserver?s_tid=ta_ans_results
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format compact;
% Define your workbook file name.
excelFullFileName = fullfile(pwd, '\testAIP.xlsx');
% Make sure it exists.  Open Excel as an ActiveX server if it does.
if isfile(excelFullFileName)
    % If the workbook exists, launch Excel as an ActiveX server.
    Excel = actxserver('Excel.Application');
    Excel.visible = true; % Make the server visible.
    fprintf('Excel opened successfully.\n');
    fprintf('Your workbook file exists:\n"%s".\nAbout to try to open it.\n', excelFullFileName);
    % Open up the existing workbook named in the variable fullFileName.
    Excel.Workbooks.Open(excelFullFileName);
    fprintf('Excel opened file successfully.\n');
else
    % File does not exist.  Alert the user.
    warningMessage = sprintf('File does not exist:\n\n"%s"\n', excelFullFileName);
    fprintf('%s\n', warningMessage);
    errordlg(warningMessage);
    return;
end
% If we get here, the workbook file exists and has been opened by Excel.
% Ask Excel for the Microsoft Azure Information Protection (AIP) label of the workbook we just opened.
label = Excel.ActiveWorkbook.SensitivityLabel.GetLabel 
% See if there is a label already.  If not, these will be null:
existingLabelID = label.LabelId 
existingLabelName = label.LabelName 
% Create a label.
label = Excel.ActiveWorkbook.SensitivityLabel.CreateLabelInfo 
label.LabelId = "a518e53f-798e-43aa-978d-c3fda1f3a682";
label.LabelName = "Business Use";
% Assign the label to the workbook.
fprintf('Setting Microsoft Azure Information Protection to "Business Use", GUID of a518e53f-798e-43aa-978d-c3fda1f3a682\n');
Excel.ActiveWorkbook.SensitivityLabel.SetLabel(label, label);
% Save this workbook with the new AIP setting we just created.
Excel.ActiveWorkbook.Save;
% Shut down Excel.
Excel.ActiveWorkbook.Close;
Excel.Quit;
% Excel is now closed down.  Delete the variable from the MATLAB workspace.
clear Excel;
% Now check to see if the AIP label has been set 
% by opening up the file in Excel and looking at the AIP banner.
winopen(excelFullFileName)

Note that there is a line in there that gets an AIP label from the existing workbook, if there is one at all. If there is not one, you can set one. But to determine what the proper LabelId (that crazy long hexadecimal number) should be, you will probably need to open an existing document that already has the label that you want set (applied to it) and then read that label with this line:

label = Excel.ActiveWorkbook.SensitivityLabel.GetLabel

John D'Errico
Last activity 2024-9-23，7:43

Using MATLAB to find a generative equation for a sequence

This stems purely from some play on my part. Suppose I asked you to work with the sequence formed as 2*n*F_n + 1, where F_n is the n'th Fibonacci number? Part of me would not be surprised to find there is nothing simple we could do. But, then it costs nothing to try, to see where MATLAB can take me in an explorative sense.

n = sym(0:100).';

Fn = fibonacci(n);

Sn = 2*n.*Fn + 1;

Sn(1:10) % A few elements

ans =

For kicks, I tried asking ChatGPT. Giving it nothing more than the first 20 members of thse sequence as integers, it decided this is a Perrin sequence, and gave me a recurrence relation, but one that is in fact incorrect. Good effort from the Ai, but a fail in the end.

Is there anything I can do? Try null! (Look carefully at the array generated by Toeplitz. It is at least a pretty way to generate the matrix I needed.)

X = toeplitz(Sn,[1,zeros(1,4)]);
rank(X(5:end,:))
ans = 5

Hmm. So there is no linear combination of those columns that yields all zeros, since the resulting matrix was full rank.

X = toeplitz(Sn,[1,zeros(1,5)]);
rank(X(6:end,:))
ans = 5

But if I take it one step further, we see the above matrix is now rank deficient. What does that tell me? It says there is some simple linear combination of the columns of X(6:end,:) that always yields zero. The previous test tells me there is no shorter constant coefficient recurrence releation, using fewer terms.

null(X(6:end,:))

ans =

Let me explain what those coefficients tell me. In fact, they yield a very nice recurrence relation for the sequence S_n, not unlike the original Fibonacci sequence it was based upon.

S(n+1) = 3*S(n) - S_(n-1) - 3*S(n-2) + S(n-3) + S(n-4)

where the first 5 members of that sequence are given as [1 3 5 13 25]. So a 6 term linear constant coefficient recurrence relation. If it reminds you of the generating relation for the Fibonacci sequence, that is good, because it should. (Remember I started the sequence at n==0, IF you decide to test it out.) We can test it out, like this:

SfunM = memoize(@(N) Sfun(N));
SfunM(25)
ans = 3751251
2*25*fibonacci(sym(25)) + 1
ans = 3751251

And indeed, it works as expected.

function Sn = Sfun(n)
switch n
    case 0
        Sn = 1;
    case 1
        Sn = 3;
    case 2
        Sn = 5;
    case 3
        Sn = 13;
    case 4
        Sn = 25;
    otherwise
        Sn = Sfun(n-5) + Sfun(n-4) - 3*Sfun(n-3) - Sfun(n-2) +3*Sfun(n-1); 
end
end

A beauty of this, is I started from nothing but a sequence of integers, derived from an expression where I had no rational expectation of finding a formula, and out drops something pretty. I might call this explorational mathematics.

The next step of course is to go in the other direction. That is, given the derived recurrence relation, if I substitute the formula for S_n in terms of the Fibonacci numbers, can I prove it is valid in general? (Yes.) After all, without some proof, it may fail for n larger than 100. (I'm not sure how much I can cram into a single discussion, so I'll stop at this point for now. If I see interest in the ideas here, I can proceed further. For example, what was I doing with that sequence in the first place? And of course, can I prove the relation is valid? Can I do so using MATLAB?)

(I'll be honest, starting from scratch, I'm not sure it would have been obvious to find that relation, so null was hugely useful here.)

Off Topic

Athanasios Paraskevopoulos
Last activity 2024-7-31

Looking for Beginner Courses on Drone Programming

Hi everyone,

I've recently joined a forest protection team in Greece, where we use drones for various tasks. This has sparked my interest in drone programming, and I'd like to learn more about it. Can anyone recommend any beginner-friendly courses or programs that teach drone programming?

I'm particularly interested in courses that focus on practical applications and might align with the work we do in forest protection. Any suggestions or guidance would be greatly appreciated!

Thank you!

lewis
Last activity 2024-7-27

"Satellite Communication simulation(main title)" Undergraduate projects, Do you have any good Ideas?

I have picked the title but don't know which direction to take it. Looking for any and all inspiration. I took the project as it sounded interesting when reading into it, but I'm a satellite novice, and my degree is in electronics.

Off Topic

Muhammad
Last activity 2024-8-26

What do you like most about MATLAB?

"What are your favorite features or functionalities in MATLAB, and how have they positively impacted your projects or research? Any tips or tricks to share?

DH
Last activity 2024-9-23，6:07

Can someone tell me what's wrong with this code? This is a problem from Cody. For example, when n=3, the output should be the matrix [1 2 2 3 3 3]

function ans = your_fcn_name(n)
n;
j=sum(1:n);
a=zeros(1,j);
 for i=1:n
   a(1,((sum(1:(i-1))+1)):(sum(1:(i-1))+i))=i.*ones(1,i);     
 end
 disp

Athanasios Paraskevopoulos
Last activity 2024-7-23

Gabriel's Horn

Gabriel's horn is a shape with the paradoxical property that it has infinite surface area, but a finite volume.

Gabriel’s horn is formed by taking the graph of

with the domain

and rotating it in three dimensions about the

axis.

There is a standard formula for calculating the volume of this shape, for a general function

.Wwe will just state that the volume of the

solid between a and b is:

The surface area of the solid is given by:

One other thing we need to consider is that we are trying to find the value of these integrals between 1 and ∞. An integral with a limit of infinity is called an improper integral and we can't evaluate it simply by plugging the value infinity into the normal equation for a definite integral. Instead, we must first calculate the definite integral up to some finite limit b and then calculate the limit of the result as b tends to ∞:

Volume

We can calculate the horn's volume using the volume integral above, so

The total volume of this infinitely long trumpet isπ.

Surface Area

To determine the surface area, we first need the function’s derivative:

Now plug it into the surface area formula and we have:

This is an improper integral and it's hard to evaluate, but since in our interval

So, we have :

Now,we evaluate this last integral

So the surface are is infinite.

% Define the function for Gabriel's Horn

gabriels_horn = @(x) 1 ./ x;

% Create a range of x values

x = linspace(1, 40, 4000); % Increase the number of points for better accuracy

y = gabriels_horn(x);

% Create the meshgrid

theta = linspace(0, 2 * pi, 6000); % Increase theta points for a smoother surface

[X, T] = meshgrid(x, theta);

Y = gabriels_horn(X) .* cos(T);

Z = gabriels_horn(X) .* sin(T);

% Plot the surface of Gabriel's Horn

figure('Position', [200, 100, 1200, 900]);

surf(X, Y, Z, 'EdgeColor', 'none', 'FaceAlpha', 0.9);

hold on;

% Plot the central axis

plot3(x, zeros(size(x)), zeros(size(x)), 'r', 'LineWidth', 2);

% Set labels

xlabel('x');

ylabel('y');

zlabel('z');

% Adjust colormap and axis properties

colormap('gray');

shading interp; % Smooth shading

% Adjust the view

view(3);

axis tight;

grid on;

% Add formulas as text annotations

dim1 = [0.4 0.7 0.3 0.2];

annotation('textbox',dim1,'String',{'$$V = \pi \int_{1}^{a} \left( \frac{1}{x} \right)^2 dx = \pi \left( 1 - \frac{1}{a} \right)$$', ...

'', ... % Add an empty line for larger gap

'$$\lim_{a \to \infty} V = \lim_{a \to \infty} \pi \left( 1 - \frac{1}{a} \right) = \pi$$'}, ...

'Interpreter','latex','FontSize',12, 'EdgeColor','none', 'FitBoxToText', 'on');

dim2 = [0.4 0.5 0.3 0.2];

annotation('textbox',dim2,'String',{'$$A = 2\pi \int_{1}^{a} \frac{1}{x} \sqrt{1 + \left( -\frac{1}{x^2} \right)^2} dx > 2\pi \int_{1}^{a} \frac{dx}{x} = 2\pi \ln(a)$$', ...

'', ... % Add an empty line for larger gap

'$$\lim_{a \to \infty} A \geq \lim_{a \to \infty} 2\pi \ln(a) = \infty$$'}, ...

'Interpreter','latex','FontSize',12, 'EdgeColor','none', 'FitBoxToText', 'on');

% Add Gabriel's Horn label

dim3 = [0.3 0.9 0.3 0.1];

annotation('textbox',dim3,'String','Gabriel''s Horn', ...

'Interpreter','latex','FontSize',14, 'EdgeColor','none', 'HorizontalAlignment', 'center');

hold off

daspect([3.5 1 1]) % daspect([x y z])

view(-27, 15)

lightangle(-50,0)

lighting('gouraud')

The properties of this figure were first studied by Italian physicist and mathematician Evangelista Torricelli in the 17th century.

Acknowledgment

I would like to express my sincere gratitude to all those who have supported and inspired me throughout this project.

First and foremost, I would like to thank the mathematician and my esteemed colleague, Stavros Tsalapatis, for inspiring me with the fascinating subject of Gabriel's Horn.

I am also deeply thankful to Mr. @Star Strider for his invaluable assistance in completing the final code.

References:

https://www.youtube.com/watch?v=z_-jSmfqlHY&t=14s

binbin
Last activity 2024-7-23

The moment the PMOS switch is turned on, the inrush current is too large and the PMOS burns out…

When it comes to MOS tube burnout, it is usually because it is not working in the SOA workspace, and there is also a case where the MOS tube is overcurrent.

For example, the maximum allowable current of the PMOS transistor in this circuit is 50A, and the maximum current reaches 80+ at the moment when the MOS transistor is turned on, then the current is very large.

At this time, the PMOS is over-specified, and we can see on the SOA curve that it is not working in the SOA range, which will cause the PMOS to be damaged.

So what if you choose a higher current PMOS? Of course you can, but the cost will be higher.

We can choose to adjust the peripheral resistance or capacitor to make the PMOS turn on more slowly, so that the current can be lowered.

For example, when adjusting R1, R2, and the jumper capacitance between gs, when Cgs is adjusted to 1uF, The Ids are only 40A max, which is fine in terms of current, and meets the 80% derating.

(50 amps * 0.8 = 40 amps).

Next, let’s look at the power, from the SOA curve, the opening time of the MOS tube is about 1ms, and the maximum power at this time is 280W.

The normal thermal resistance of the chip is 50°C/W, and the maximum junction temperature can be 302°F.

Assuming the ambient temperature is 77°F, then the instantaneous power that 1ms can withstand is about 357W.

The actual power of PMOS here is 280W, which does not exceed the limit, which means that it works normally in the SOA area.

Therefore, when the current impact of the MOS transistor is large at the moment of turning, the Cgs capacitance can be adjusted appropriately to make the PMOS Working in the SOA area, you can avoid the problem of MOS corruption.

Marisa
Last activity 2024-8-26

Can't click Roll the Dice!

I am trying to earn my Intro to MATLAB badge in Cody, but I cannot click the Roll the Dice! problem. It simply is not letting me click it, therefore I cannot earn my badge. Does anyone know who I should contact or what to do?

Muhammad
Last activity 2024-7-27

Need help to get started with matlab

Hello everyone, i hope you all are in good health. i need to ask you about the help about where i should start to get indulge in matlab. I am an electrical engineer but having experience of construction field. I am new here. Please do help me. I shall be waiting forward to hear from you. I shall be grateful to you. Need recommendations and suggestions from experienced members. Thank you.

James Blanchard
Last activity 2024-7-2

document on solving ODEs and PDEs

I recently wrote up a document which addresses the solution of ordinary and partial differential equations in Matlab (with some Python examples thrown in for those who are interested). For ODEs, both initial and boundary value problems are addressed. For PDEs, it addresses parabolic and elliptic equations. The emphasis is on finite difference approaches and built-in functions are discussed when available. Theory is kept to a minimum. I also provide a discussion of strategies for checking the results, because I think many students are too quick to trust their solutions. For anyone interested, the document can be found at https://blanchard.neep.wisc.edu/SolvingDifferentialEquationsWithMatlab.pdf

Sylvester
Last activity 2024-6-26

Geometry-Based Stochastic Channel Modeling

Kindly link me to the Channel Modeling Group.

I read and compreheneded a paper on channel modeling "An Adaptive Geometry-Based Stochastic Model for Non-Isotropic MIMO Mobile-to-Mobile Channels" except the graphical results obtained from the MATLAB codes. I have tried to replicate the same graphs but to no avail from my codes. And I am really interested in the topic, i have even written to the authors of the paper but as usual, there is no reply from them. Kindly assist if possible.

wireless channel modeling

Saleem
Last activity 2024-6-27

Where to find coding and algorithms problems with solutions?

Hi, I'm looking for sites where I can find coding & algorithms problems and their solutions. I'm doing this workshop in college and I'll need some problems to go over with the students and explain how Matlab works by solving the problems with them and then reviewing and going over different solution options. Does anyone know a website like that? I've tried looking in the Matlab Cody By Mathworks, but didn't exactly find what I'm looking for. Thanks in advance.

Athanasios Paraskevopoulos
Last activity 2024-7-4

Stochastic simulation of a novel pathogen

We are modeling the introduction of a novel pathogen into a completely susceptible population. In the cells below, I have provided you with the Matlab code for a simple stochastic SIR model, implemented using the "GillespieSSA" function

Simulating the stochastic model 100 times for

Since γ is 0.4 per day,

per day

% Define the parameters

beta = 0.36;

gamma = 0.4;

n_sims = 100;

tf = 100; % Time frame changed to 100

% Calculate R0

R0 = beta / gamma

R0 = 0.9000

% Initial state values

initial_state_values = [1000000; 1; 0; 0]; % S, I, R, cum_inc

% Define the propensities and state change matrix

a = @(state) [beta * state(1) * state(2) / 1000000, gamma * state(2)];

nu = [-1, 0; 1, -1; 0, 1; 0, 0];

% Define the Gillespie algorithm function

function [t_values, state_values] = gillespie_ssa(initial_state, a, nu, tf)

t = 0;

state = initial_state(:); % Ensure state is a column vector

t_values = t;

state_values = state';

while t < tf

rates = a(state);

rate_sum = sum(rates);

if rate_sum == 0

break;

end

tau = -log(rand) / rate_sum;

t = t + tau;

r = rand * rate_sum;

cum_sum_rates = cumsum(rates);

reaction_index = find(cum_sum_rates >= r, 1);

state = state + nu(:, reaction_index);

% Update cumulative incidence if infection occurred

if reaction_index == 1

state(4) = state(4) + 1; % Increment cumulative incidence

end

t_values = [t_values; t];

state_values = [state_values; state'];

end

% Function to simulate the stochastic model multiple times and plot results

function simulate_stoch_model(beta, gamma, n_sims, tf, initial_state_values, R0, plot_type)

% Define the propensities and state change matrix

a = @(state) [beta * state(1) * state(2) / 1000000, gamma * state(2)];

nu = [-1, 0; 1, -1; 0, 1; 0, 0];

% Set random seed for reproducibility

rng(11);

% Initialize plot

figure;

hold on;

for i = 1:n_sims

[t, output] = gillespie_ssa(initial_state_values, a, nu, tf);

% Check if the simulation had only one step and re-run if necessary

while length(t) == 1

[t, output] = gillespie_ssa(initial_state_values, a, nu, tf);

end

if strcmp(plot_type, 'cumulative_incidence')

plot(t, output(:, 4), 'LineWidth', 2, 'Color', rand(1, 3));

elseif strcmp(plot_type, 'prevalence')

plot(t, output(:, 2), 'LineWidth', 2, 'Color', rand(1, 3));

end

xlabel('Time (days)');

if strcmp(plot_type, 'cumulative_incidence')

ylabel('Cumulative Incidence');

ylim([0 inf]);

elseif strcmp(plot_type, 'prevalence')

ylabel('Prevalence of Infection');

ylim([0 50]);

end

title(['Stochastic model output for R0 = ', num2str(R0)]);

subtitle([num2str(n_sims), ' simulations']);

xlim([0 tf]);

grid on;

hold off;

end

% Simulate the model 100 times and plot cumulative incidence

simulate_stoch_model(beta, gamma, n_sims, tf, initial_state_values, R0, 'cumulative_incidence');

% Simulate the model 100 times and plot prevalence

simulate_stoch_model(beta, gamma, n_sims, tf, initial_state_values, R0, 'prevalence');

https://blogs.mathworks.com/community/2015/02/26/lets-code-make-a-cody-video/

Ned Gulley
Last activity 2024-6-13

Make a video: "How I Solved It"

Twitch built an entire business around letting you watch over someone's shoulder while they play video games. I feel like we should be able to make at least a few videos where we get to watch over someone's shoulder while they solve Cody problems. I would pay good money for a front-row seat to watch some of my favorite solvers at work. Like, I want to know, did Alfonso Nieto-Castonon just sit down and bang out some of those answers, or did he have to think about it for a while? What was he thinking about while he solved it? What resources was he drawing on? There's nothing like watching a master craftsman at work.

I can imagine a whole category of Cody videos called "How I Solved It". I tried making one of these myself a while back, but as far as I could tell, nobody else made one.

Here's the direct link to the video: https://www.youtube.com/watch?v=hoSmO1XklAQ

I hereby challenge you to make a "How I Solved It" video and post it here. If you make one, I'll make another one.

goc3
Last activity 2024-9-8，14:47

+1 for Backwards Loops!

Base case:

Suppose you need to do a computation many times. We are going to assume that this computation cannot be vectorized. The simplest case is to use a for loop:

number_of_elements = 1e6;
test_fcn = @(x) sqrt(x) / x;
tic
for i = 1:number_of_elements
    x(i) = test_fcn(i);
end
t_forward = toc;
disp(t_forward + " seconds")
0.10925 seconds

Preallocation:

This can easily be sped up by preallocating the variable that houses results:

tic
x = zeros(number_of_elements, 1);
for i = 1:number_of_elements
    x(i) = test_fcn(i);
end
t_forward_prealloc = toc;
disp(t_forward_prealloc + " seconds")
0.035106 seconds

In this example, preallocation speeds up the loop by a factor of about three to four (running in R2024a). Comment below if you get dramatically different results.

disp(sprintf("%.1f", t_forward / t_forward_prealloc))
3.1

Run it in reverse:

Is there a way to skip the explicit preallocation and still be fast? Indeed, there is.

clear x
tic
for i = number_of_elements:-1:1
    x(i) = test_fcn(i);
end
t_backward = toc;
disp(t_backward + " seconds")
0.032392 seconds

By running the loop backwards, the preallocation is implicitly performed during the first iteration and the loop runs in about the same time (within statistical noise):

disp(sprintf("%.2f", t_forward_prealloc / t_backward))
1.08

Do you get similar results when running this code? Let us know your thoughts in the comments below.

Beneficial side effect:

Have you ever had to use a for loop to delete elements from a vector? If so, keeping track of index offsets can be tricky, as deleting any element shifts all those that come after. By running the for loop in reverse, you don't need to worry about index offsets while deleting elements.