Search Discussions

Athanasios Paraskevopoulos
Last activity 2024 年 8 月 12 日

Strange attractors

An attractor is called strange if it has a fractal structure, that is if it has non-integer Hausdorff dimension. This is often the case when the dynamics on it are chaotic, but strange nonchaotic attractors also exist. If a strange attractor is chaotic, exhibiting sensitive dependence on initial conditions, then any two arbitrarily close alternative initial points on the attractor, after any of various numbers of iterations, will lead to points that are arbitrarily far apart (subject to the confines of the attractor), and after any of various other numbers of iterations will lead to points that are arbitrarily close together. Thus a dynamic system with a chaotic attractor is locally unstable yet globally stable: once some sequences have entered the attractor, nearby points diverge from one another but never depart from the attractor.

The term strange attractor was coined by David Ruelle and Floris Takens to describe the attractor resulting from a series of bifurcations of a system describing fluid flow. Strange attractors are often differentiable in a few directions, but some are like a Cantor dust, and therefore not differentiable. Strange attractors may also be found in the presence of noise, where they may be shown to support invariant random probability measures of Sinai–Ruelle–Bowen type.

Examples of strange attractors include the Rössler attractor, and Lorenz attractor.

Lorenz

% Lorenz Attractor Parameters

sigma = 10;

beta = 8/3;

rho = 28;

% Lorenz system of differential equations

f = @(t, a) [-sigma*a(1) + sigma*a(2);

rho*a(1) - a(2) - a(1)*a(3);

-beta*a(3) + a(1)*a(2)];

% Time span

tspan = [0 100];

% Initial conditions

a0 = [1 1 1];

% Solve the system using ode45

[t, a] = ode45(f, tspan, a0);

% Plot using scatter3 with time-based color mapping

figure;

scatter3(a(:,1), a(:,2), a(:,3), 5, t, 'filled'); % 5 is the marker size

title('Lorenz Attractor');

xlabel('x(t)');

ylabel('y(t)');

zlabel('z(t)');

grid on;

colorbar; % Add a colorbar to indicate the time mapping

view(3); % Set the view to 3D

Sprott

% Define the parameters

a = 2.07;

b = 1.79;

% Define the system of differential equations

dynamics = @(t, X) [ ...

X(2) + a * X(1) * X(2) + X(1) * X(3); % dx/dt

1 - b * X(1)^2 + X(2) * X(3); % dy/dt

X(1) - X(1)^2 - X(2)^2 % dz/dt

];

% Initial conditions

X0 = [0.63; 0.47; -0.54];

% Time span

tspan = [0 100];

% Solve the system using ode45

[t, X] = ode45(dynamics, tspan, X0);

% Plot the results with color gradient

figure;

colormap(jet); % Set the colormap

c = linspace(1, 10, length(t)); % Color data based on time

% Create a 3D line plot with color based on time

for i = 1:length(t)-1

plot3(X(i:i+1,1), X(i:i+1,2), X(i:i+1,3), 'Color', [0 0.5 0.9]*c(i)/10, 'LineWidth', 1.5);

hold on;

end

% Set plot properties

title('Sprott Attractor');

xlabel('x(t)');

ylabel('y(t)');

zlabel('z(t)');

grid on;

colorbar; % Add a colorbar to indicate the time mapping

view(3); % Set the view to 3D

hold off;

Rössler

% Define the parameters

a = 0.2;

b = 0.2;

c = 5.7;

% Define the system of differential equations

dynamics = @(t, X) [ ...

-(X(2) + X(3)); % dx/dt

X(1) + a * X(2); % dy/dt

b + X(3) * (X(1) - c) % dz/dt

];

% Initial conditions

X0 = [10.0; 0.00; 10.0];

% Time span

tspan = [0 100];

% Solve the system using ode45

[t, X] = ode45(dynamics, tspan, X0);

% Plot the results

figure;

scatter3(X(:,1), X(:,2), X(:,3), 5, t, 'filled');

title('Rössler Attractor');

xlabel('x(t)');

ylabel('y(t)');

zlabel('z(t)');

grid on;

colorbar; % Add a colorbar to indicate the time mapping

view(3); % Set the view to 3D

Rabinovich-Fabrikant

%% Parameters for Rabinovich-Fabrikant Attractor

alpha = 0.14;

gamma = 0.10;

dt = 0.01;

num_steps = 5000;

% Initial conditions

x0 = -1;

y0 = 0;

z0 = 0.5;

% Preallocate arrays for performance

x = zeros(1, num_steps);

y = zeros(1, num_steps);

z = zeros(1, num_steps);

% Set initial values

x(1) = x0;

y(1) = y0;

z(1) = z0;

% Generate the attractor

for i = 1:num_steps-1

x(i+1) = x(i) + dt * (y(i)*(z(i) - 1 + x(i)^2) + gamma*x(i));

y(i+1) = y(i) + dt * (x(i)*(3*z(i) + 1 - x(i)^2) + gamma*y(i));

z(i+1) = z(i) + dt * (-2*z(i)*(alpha + x(i)*y(i)));

end

% Create a time vector for color mapping

t = linspace(0, 100, num_steps);

% Plot using scatter3

figure;

scatter3(x, y, z, 5, t, 'filled'); % 5 is the marker size

title('Rabinovich-Fabrikant Attractor');

xlabel('x(t)');

ylabel('y(t)');

zlabel('z(t)');

grid on;

colorbar; % Add a colorbar to indicate the time mapping

view(3); % Set the view to 3D

References

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 年 9 月 9 日

What Olympic sports to you like to watch most?

Swimming, diving

16%

Other water-based sport

Gymnastics

20%

Other indoor arena sport

15%

track, field

24%

Other outdoor sport

21%

346 票

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 日

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.)

How can i edit my code which works on r2014b version at work but not on my personal r2024a version? by Oluwadamilola Oke

David
Last activity 2024 年 7 月 30 日

Week of July 30th - Must-See MATLAB Central Posts

Hello, MATLAB enthusiasts! 🌟

Over the past few weeks, our community has been buzzing with insightful questions, vibrant discussions, and innovative ideas. Whether you're a seasoned expert or a curious beginner, there's something here for everyone to learn and enjoy. Let's take a moment to highlight some of the standout contributions that have sparked interest and inspired many. Dive in and see how you can join the conversation or find solutions to your own challenges!

Interesting Questions

Oluwadamilola Oke is seeking assistance with a MATLAB code that works on version r2014b but encounters errors on version r2024a. The issue seems to be related to file location or the use of specific commands like movefile. If you have experience with these versions of MATLAB, your expertise could be invaluable.

Something doesn't work for me in fit by Yohay

Yohay has been working on a simulation to measure particle speed and fit it to the Maxwell-Boltzmann distribution. However, the fit isn't aligning perfectly with the data. Yohay has shared the code and histogram data for community members to review and provide suggestions.

Request about comment char % vs C++ // by Alessandro Livi

Alessandro Livi is toggling between C++ for Arduino Pico and MATLAB App Designer. They suggest an enhancement where typing // for comments in MATLAB automatically converts to %. This small feature could improve the workflow for many users who switch between programming languages.

Popular Discussions

Gabriel's Horn by Athanasios Paraskevopoulos

Athanasios Paraskevopoulos has started an engaging discussion on Gabriel's Horn, a shape with infinite surface area but finite volume. The conversation delves into the mathematical intricacies and integral calculations required to understand this paradoxical shape.

Matlab Custom Font by Honzik

Honzik has brought up an interesting topic about custom fonts for MATLAB. While popular coding fonts handle characters like 0 and O well, they often fail to distinguish between different types of brackets. Honzik suggests that MathWorks could develop a custom font optimized for MATLAB syntax to reduce coding errors.

From the Blogs

Derivative of state ‘1’ in block ‘X/Y/Integrator’ at time 0.55 is not finite By Guy Rouleau

Guy Rouleau addresses a common error in Simulink models: "Derivative of state '1' in block 'X/Y/Integrator' at time 0.55 is not finite." The blog post explores various tools and methods to diagnose and resolve this issue, making it a valuable read for anyone facing similar challenges.

Animating Science: Creating Time-Sensitive Animations with MATLAB by Gianluca Carnielli

Guest writer Gianluca Carnielli, featured by Adam Danz, shares insights on creating time-sensitive animations using MATLAB. The article covers controlling the motion of multiple animated objects, organizing data with timetables, and simplifying animations with the retime function. This is a must-read for anyone interested in scientific animations.

Feel free to check out these fascinating contributions and join the discussions! Your input and expertise can make a significant difference in our community.

must-see

supercomputers
Last activity 2024 年 11 月 6 日 1:25

writing matlab code using llm

hello i found the following tools helpful to write matlab programs. copilot.microsoft.com chatgpt.com/gpts gemini.google.com and ai.meta.com. thanks a lot and best wishes.

https://www.mathworks.com/matlabcentral/fileexchange/163796-large-language-models-llms-with-matlab

Hans Scharler
Last activity 2024 年 7 月 26 日

LLMs with MATLAB Now Supports OpenAI GPT-4o mini

Check out the LLMs with MATLAB project on File Exchange to access Large Language Models from MATLAB.

Along with the latest support for GPT-4o mini, you can use LLMs with MATLAB to generate images, categorize data, and provide semantic analyis.

Run it now by clicking Open in MATLAB Online, signing in, and using your API Key from OpenAI.

Kalhara
Last activity 2024 年 6 月 19 日

NVIDIA supercharge the AI world

What do you think about the NVIDIA's achivement of becoming the top giant of manufacturing chips, especially for AI world?

Image Analyst
Last activity 2024 年 11 月 6 日 2:26

Ask Me Anything about image analysis or the Mathworks community

Hello, everyone! I’m Mark Hayworth, but you might know me better in the community as Image Analyst. I've been using MATLAB since 2006 (18 years). My background spans a rich career as a former senior scientist and inventor at The Procter & Gamble Company (HQ in Cincinnati). I hold both master’s & Ph.D. degrees in optical sciences from the College of Optical Sciences at the University of Arizona, specializing in imaging, image processing, and image analysis. I have 40+ years of military, academic, and industrial experience with image analysis programming and algorithm development. I have experience designing custom light booths and other imaging systems. I also work with color and monochrome imaging, video analysis, thermal, ultraviolet, hyperspectral, CT, MRI, radiography, profilometry, microscopy, NIR, and Raman spectroscopy, etc. on a huge variety of subjects.

I'm thrilled to participate in MATLAB Central's Ask Me Anything (AMA) session, a fantastic platform for knowledge sharing and community engagement. Following Adam Danz’s insightful AMA on staff contributors in the Answers forum, I’d like to discuss topics in the area of image analysis and processing. I invite you to ask me anything related to this field, whether you're seeking recommendations on tools, looking for tips and tricks, my background, or career development advice. Additionally, I'm more than willing to share insights from my experiences in the MATLAB Answers community, File Exchange, and my role as a member of the Community Advisory Board. If you have questions related to your specific images or your custom MATLAB code though, I'll invite you to ask those in the Answers forum. It's a more appropriate forum for those kinds of questions, plus you can get the benefit of other experts offering their solutions in addition to me.

For the coming weeks, I'll be here to engage with your questions and help shed light on any topics you're curious about.

108

Proving one function is greater than other? by Fatima Majeed

David
Last activity 2024 年 6 月 17 日

Week of June 17th - Must-See MATLAB Central Posts

Hello, everyone!

Over the past few weeks, our community has been buzzing with activity, showcasing the incredible depth of knowledge, creativity, and innovation that makes this forum such a vibrant place. Today, we're excited to highlight some of the noteworthy contributions that have sparked discussions, offered insights, and shared knowledge across various topics. Let's dive in!

Interesting Questions

Fatima Majeed brings us a thought-provoking mathematical challenge, delving into inequalities and the realms beyond (e^e). If you're up for a mathematical journey, this question is a must-see!

How to split/segment csv file into multiple sections? by lil brain

lil brain tackles a practical problem many of us have faced: efficiently segmenting a CSV file based on specific criteria. This post is not only a query but a learning opportunity for anyone dealing with similar data manipulation challenges.

Popular Discussions

Trick: Easy Digit Manipulation by goc3

Discover a simple yet effective trick for digit manipulation from goc3. This tip is especially handy for those frequenting Cody challenges or anyone interested in enhancing their number handling skills in MATLAB.

Run Code in Discussions Area! by Chen Lin

Chen Lin shares an exciting update about the 'Run Code' feature in the Discussions area, highlighting how our community can now directly execute and share code snippets within discussions. This feature marks a significant enhancement in how we interact and solve problems together.

From the Blogs

A Deep Dive into EEG Analysis for Predicting Neurological Outcomes By Tanya Kuruvilla

Connell D`Souza, alongside Team Swarthbeat, explores the cutting-edge application of EEG analysis in predicting neurological outcomes post-cardiac arrest. This blog post offers an in-depth look into the challenges and methodologies of modern medical data analysis.

Crafting the Robots of Tomorrow: The Power of Robot Simulation with MATLAB and Simulink by YJ Lim

Mihir Acharya discusses the pivotal role of MATLAB and Simulink in the future of robotics simulation. Through an engaging conversation with industry analyst George Chowdhury, this post sheds light on overcoming simulation challenges and the exciting possibilities that lie ahead.

We encourage everyone to explore these contributions further and engage with the authors and the community. Your participation is what fuels this community's continual growth and innovation.

Here's to many more discussions, discoveries, and breakthroughs together!

must-see

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');

goc3
Last activity 2024 年 9 月 8 日

+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.

Chen Lin
Last activity 2024 年 6 月 14 日

Run Code in Discussions Area!

We're thrilled to share an exciting update with our community: the 'Run Code' feature is now available in the Discussions area!

Simply insert your code into the editor and press the green triangle button to run it. Your code will execute using the latest MATLAB R24a version, and it supports most common toolboxes. Moreover, this innovative feature allows for the running of attached files, further enhancing its utility and flexibility.

The ‘run code’ feature was first introduced in MATLAB Answers. Encouraged by the positive feedback and at the request of our community members, we are now expanding the availability of this feature to more areas within our community.

As always, your feedback is crucial to us, so please don't hesitate to share your thoughts and experiences by leaving a comment.

Zhaoxu Liu / slandarer
Last activity 2024 年 6 月 9 日

How to set the AlphaData of a colorbar?

Many times when ploting, we not only need to set the color of the plot, but also its

transparency, Then how we set the alphaData of colorbar at the same time ?

It seems easy to do so :

data = rand(12,12); 
% Transparency range 0-1, .3-1 for better appearance here
AData = rescale(- data, .3, 1);
% Draw an imagesc with numerical control over colormap and transparency
imagesc(data, 'AlphaData',AData);
colormap(jet); 
ax = gca;
ax.DataAspectRatio = [1,1,1];
ax.TickDir = 'out';          
ax.Box = 'off';
% get colorbar object
CBarHdl = colorbar;
pause(1e-16)
% Modify the transparency of the colorbar
CData = CBarHdl.Face.Texture.CData;
ALim = [min(min(AData)), max(max(AData))];
CData(4,:) = uint8(255.*rescale(1:size(CData, 2), ALim(1), ALim(2)));
CBarHdl.Face.Texture.ColorType = 'TrueColorAlpha';
CBarHdl.Face.Texture.CData = CData;

But !!!!!!!!!!!!!!! We cannot preserve the changes when saving them as images ：

It seems that when saving plots, the `Texture` will be refresh, but the `Face` will not :

however, object Face only have 4 colors to change(The four corners of a quadrilateral), how

can we set more colors ??

`Face` is a quadrilateral object, and we can change the `VertexData` to draw more than one little quadrilaterals:

data = rand(12,12); 
% Transparency range 0-1, .3-1 for better appearance here
AData = rescale(- data, .3, 1);
%Draw an imagesc with numerical control over colormap and transparency
imagesc(data, 'AlphaData',AData);
colormap(jet); 
 
ax = gca;
ax.DataAspectRatio = [1,1,1];  
ax.TickDir = 'out';           
ax.Box = 'off';
% get colorbar object
CBarHdl = colorbar;
pause(1e-16)
% Modify the transparency of the colorbar
CData = CBarHdl.Face.Texture.CData;
ALim = [min(min(AData)), max(max(AData))];
CData(4,:) = uint8(255.*rescale(1:size(CData, 2), ALim(1), ALim(2)));
warning off
CBarHdl.Face.ColorType = 'TrueColorAlpha';
VertexData = CBarHdl.Face.VertexData;
tY = repmat((1:size(CData,2))./size(CData,2), [4,1]);
tY1 = tY(:).'; tY2 = tY  - tY(1,1); tY2(3:4,:) = 0; tY2 = tY2(:).';
tM1 = [tY1.*0 + 1; tY1; tY1.*0 + 1];
tM2 = [tY1.*0; tY2; tY1.*0];
CBarHdl.Face.VertexData = repmat(VertexData, [1,size(CData,2)]).*tM1 + tM2;
CBarHdl.Face.ColorData = reshape(repmat(CData, [4,1]), 4, []);

The higher the value, the more transparent it becomes

data = rand(12,12); 
AData = rescale(- data, .3, 1);
imagesc(data, 'AlphaData',AData);
colormap(jet); 
ax = gca;
ax.DataAspectRatio = [1,1,1];
ax.TickDir = 'out';          
ax.Box = 'off';
CBarHdl = colorbar;
pause(1e-16)
CData = CBarHdl.Face.Texture.CData;
ALim = [min(min(AData)), max(max(AData))];
CData(4,:) = uint8(255.*rescale(size(CData, 2):-1:1, ALim(1), ALim(2)));
warning off
CBarHdl.Face.ColorType = 'TrueColorAlpha';
VertexData = CBarHdl.Face.VertexData;
tY = repmat((1:size(CData,2))./size(CData,2), [4,1]);
tY1 = tY(:).'; tY2 = tY  - tY(1,1); tY2(3:4,:) = 0; tY2 = tY2(:).';
tM1 = [tY1.*0 + 1; tY1; tY1.*0 + 1];
tM2 = [tY1.*0; tY2; tY1.*0];
CBarHdl.Face.VertexData = repmat(VertexData, [1,size(CData,2)]).*tM1 + tM2;
CBarHdl.Face.ColorData = reshape(repmat(CData, [4,1]), 4, []);

More transparent in the middle

data = rand(12,12) - .5; 
AData = rescale(abs(data), .1, .9);
imagesc(data, 'AlphaData',AData);
colormap(jet); 
ax = gca;
ax.DataAspectRatio = [1,1,1];
ax.TickDir = 'out';          
ax.Box = 'off';
CBarHdl = colorbar;
pause(1e-16)
CData = CBarHdl.Face.Texture.CData;
ALim = [min(min(AData)), max(max(AData))];
CData(4,:) = uint8(255.*rescale(abs((1:size(CData, 2)) - (1 + size(CData, 2))/2), ALim(1), ALim(2)));
warning off
CBarHdl.Face.ColorType = 'TrueColorAlpha';
VertexData = CBarHdl.Face.VertexData;
tY = repmat((1:size(CData,2))./size(CData,2), [4,1]);
tY1 = tY(:).'; tY2 = tY  - tY(1,1); tY2(3:4,:) = 0; tY2 = tY2(:).';
tM1 = [tY1.*0 + 1; tY1; tY1.*0 + 1];
tM2 = [tY1.*0; tY2; tY1.*0];
CBarHdl.Face.VertexData = repmat(VertexData, [1,size(CData,2)]).*tM1 + tM2;
CBarHdl.Face.ColorData = reshape(repmat(CData, [4,1]), 4, []);

The code will work if the plot have AlphaData property

data = peaks(30);
AData = rescale(data, .2, 1);
surface(data, 'FaceAlpha','flat','AlphaData',AData);
colormap(jet(100)); 
ax = gca;
ax.DataAspectRatio = [1,1,1];
ax.TickDir = 'out';          
ax.Box = 'off';
view(3)
CBarHdl = colorbar;
pause(1e-16)
CData = CBarHdl.Face.Texture.CData;
ALim = [min(min(AData)), max(max(AData))];
CData(4,:) = uint8(255.*rescale(1:size(CData, 2), ALim(1), ALim(2)));
warning off
CBarHdl.Face.ColorType = 'TrueColorAlpha';
VertexData = CBarHdl.Face.VertexData;
tY = repmat((1:size(CData,2))./size(CData,2), [4,1]);
tY1 = tY(:).'; tY2 = tY  - tY(1,1); tY2(3:4,:) = 0; tY2 = tY2(:).';
tM1 = [tY1.*0 + 1; tY1; tY1.*0 + 1];
tM2 = [tY1.*0; tY2; tY1.*0];
CBarHdl.Face.VertexData = repmat(VertexData, [1,size(CData,2)]).*tM1 + tM2;
CBarHdl.Face.ColorData = reshape(repmat(CData, [4,1]), 4, []);

Books

Athanasios Paraskevopoulos
Last activity 2024 年 6 月 9 日

Discovering an Excellent Resource on Ordinary Differential Equations

While searching the internet for some books on ordinary differential equations, I came across a link that I believe is very useful for all math students and not only. If you are interested in ODEs, it's worth taking the time to study it.

A First Look at Ordinary Differential Equations by Timothy S. Judson is an excellent resource for anyone looking to understand ODEs better. Here's a brief overview of the main topics covered:

Introduction to ODEs: Basic concepts, definitions, and initial differential equations.
Methods of Solution:

Separable equations
First-order linear equations
Exact equations
Transcendental functions

Applications of ODEs: Practical examples and applications in various scientific fields.
Systems of ODEs: Analysis and solutions of systems of differential equations.
Series and Numerical Methods: Use of series and numerical methods for solving ODEs.

This book provides a clear and comprehensive introduction to ODEs, making it suitable for students and new researchers in mathematics. If you're interested, you can explore the book in more detail here: A First Look at Ordinary Differential Equations.

Chen Lin
Last activity 2024 年 6 月 13 日

Introducing a New Discussions Area for Cody Users!

Many MATLAB enthusiasts come Cody to sharpen their skills, face new challenges, and engage in friendly competition. We firmly believe that learning from peers is one of the most effective ways to grow.

With this in mind, the Cody team is thrilled to unveil a new feature aimed at enriching your learning journey: the Cody Discussion Channel. This space is designed for sharing expertise, acquiring new skills, and fostering connections within our community.

On the Cody homepage, you'll now notice a Discussions section, prominently displaying the four most recent posts. For those eager to contribute, we encourage you to familiarize yourself with our posting guidelines before creating a new post. This will help maintain a constructive and valuable exchange of ideas for everyone involved.

Together, let's create an environment where every member feels empowered to share, learn, and connect.

cody

Hans Scharler
Last activity 2024 年 5 月 30 日

Surprised that I was able to create a dice rolling simulation

One of the starter prompts is about rolling two six-sided dice and plot the results. As a hobby, I create my own board games. I was able to use the dice rolling prompt to show how a simple roll and move game would work. That was a great surprise!

playground