ans =
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
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,:))
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,:))
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,:))
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)
2*25*fibonacci(sym(25)) + 1
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.)
Following on from my previous post The Non-Chaotic Duffing Equation, now we will study the chaotic behaviour of the Duffing Equation
P.s:Any comments/advice on improving the code is welcome.
The Original Duffing Equation is the following:
Let 
. This implies that
. This implies thatThen we rewrite it as a System of First-Order Equations
Using the substitution 
for 
, the second-order equation can be transformed into the following system of first-order equations:
for , the second-order equation can be transformed into the following system of first-order equations:
Exploring the Effect of γ.
% Define parameters
gamma = 0.1;
alpha = -1;
beta = 1;
delta = 0.1;
omega = 1.4;
% Define the system of equations
odeSystem = @(t, y) [y(2);
-delta*y(2) - alpha*y(1) - beta*y(1)^3 + gamma*cos(omega*t)];
% Initial conditions
y0 = [0; 0]; % x(0) = 0, v(0) = 0
% Time span
tspan = [0 200];
% Solve the system
[t, y] = ode45(odeSystem, tspan, y0);
% Plot the results
figure;
plot(t, y(:, 1));
xlabel('Time');
ylabel('x(t)');
title('Solution of the nonlinear system');
grid on;
% Plot the phase portrait
figure;
plot(y(:, 1), y(:, 2));
xlabel('x(t)');
ylabel('v(t)');
title('Phase Portrait');
grid on;
% Define the tail (e.g., last 10% of the time interval)
tail_start = floor(0.9 * length(t)); % Starting index for the tail
tail_end = length(t); % Ending index for the tail
% Plot the tail of the solution
figure;
plot(y(tail_start:tail_end, 1), y(tail_start:tail_end, 2), 'r', 'LineWidth', 1.5);
xlabel('x(t)');
ylabel('v(t)');
title('Phase Portrait - Tail of the Solution');
grid on;
% Define parameters
gamma = 0.318;
alpha = -1;
beta = 1;
delta = 0.1;
omega = 1.4;
% Define the system of equations
odeSystem = @(t, y) [y(2);
-delta*y(2) - alpha*y(1) - beta*y(1)^3 + gamma*cos(omega*t)];
% Initial conditions
y0 = [0; 0]; % x(0) = 0, v(0) = 0
% Time span
tspan = [0 800];
% Solve the system
[t, y] = ode45(odeSystem, tspan, y0);
% Plot the results
figure;
plot(t, y(:, 1));
xlabel('Time');
ylabel('x(t)');
title('Solution of the nonlinear system');
grid on;
% Plot the phase portrait
figure;
plot(y(:, 1), y(:, 2));
xlabel('x(t)');
ylabel('v(t)');
title('Phase Portrait');
grid on;
% Define the tail (e.g., last 10% of the time interval)
tail_start = floor(0.9 * length(t)); % Starting index for the tail
tail_end = length(t); % Ending index for the tail
% Plot the tail of the solution
figure;
plot(y(tail_start:tail_end, 1), y(tail_start:tail_end, 2), 'r', 'LineWidth', 1.5);
xlabel('x(t)');
ylabel('v(t)');
title('Phase Portrait - Tail of the Solution');
grid on;
% Define parameters
gamma = 0.338;
alpha = -1;
beta = 1;
delta = 0.1;
omega = 1.4;
% Define the system of equations
odeSystem = @(t, y) [y(2);
-delta*y(2) - alpha*y(1) - beta*y(1)^3 + gamma*cos(omega*t)];
% Initial conditions
y0 = [0; 0]; % x(0) = 0, v(0) = 0
% Time span with more points for better resolution
tspan = linspace(0, 200,2000); % Increase the number of points
% Solve the system
[t, y] = ode45(odeSystem, tspan, y0);
% Plot the results
figure;
plot(t, y(:, 1));
xlabel('Time');
ylabel('x(t)');
title('Solution of the nonlinear system');
grid on;
% Plot the phase portrait
figure;
plot(y(:, 1), y(:, 2));
xlabel('x(t)');
ylabel('v(t)');
title('Phase Portrait');
grid on;
% Define the tail (e.g., last 10% of the time interval)
tail_start = floor(0.9 * length(t)); % Starting index for the tail
tail_end = length(t); % Ending index for the tail
% Plot the tail of the solution
figure;
plot(y(tail_start:tail_end, 1), y(tail_start:tail_end, 2), 'r', 'LineWidth', 1.5);
xlabel('x(t)');
ylabel('v(t)');
title('Phase Portrait - Tail of the Solution');
grid on;
ax = gca;
chart = ax.Children(1);
datatip(chart,0.5581,-0.1126);
% Define parameters
gamma = 0.35;
alpha = -1;
beta = 1;
delta = 0.1;
omega = 1.4;
% Define the system of equations
odeSystem = @(t, y) [y(2);
-delta*y(2) - alpha*y(1) - beta*y(1)^3 + gamma*cos(omega*t)];
% Initial conditions
y0 = [0; 0]; % x(0) = 0, v(0) = 0
% Time span with more points for better resolution
tspan = linspace(0, 400,3000); % Increase the number of points
% Solve the system
[t, y] = ode45(odeSystem, tspan, y0);
% Plot the results
figure;
plot(t, y(:, 1));
xlabel('Time');
ylabel('x(t)');
title('Solution of the nonlinear system');
grid on;
% Plot the phase portrait
figure;
plot(y(:, 1), y(:, 2));
xlabel('x(t)');
ylabel('v(t)');
title('Phase Portrait');
grid on;
% Define the tail (e.g., last 10% of the time interval)
tail_start = floor(0.9 * length(t)); % Starting index for the tail
tail_end = length(t); % Ending index for the tail
% Plot the tail of the solution
figure;
plot(y(tail_start:tail_end, 1), y(tail_start:tail_end, 2), 'r', 'LineWidth', 1.5);
xlabel('x(t)');
ylabel('v(t)');
title('Phase Portrait - Tail of the Solution');
grid on;
Studying the attached document Duffing Equation from the University of Colorado, I noticed that there is an analysis of The Non-Chaotic Duffing Equation and all the graphs were created with Matlab. And since the code is not given I took the initiative to try to create the same graphs with the following code.
- Plotting the Potential Energy and Identifying Extrema
% Define the range of x values
x = linspace(-2, 2, 1000);
% Define the potential function V(x)
V = -x.^2 / 2 + x.^4 / 4;
% Plot the potential function
figure;
plot(x, V, 'LineWidth', 2);
hold on;
% Mark the minima at x = ±1
plot([-1, 1], [-1/4, -1/4], 'ro', 'MarkerSize', 5, 'MarkerFaceColor', 'g');
% Add LaTeX title and labels
title('Duffing Potential Energy: $$V(x) = -\frac{x^2}{2} + \frac{x^4}{4}$$', 'Interpreter', 'latex');
xlabel('$$x$$', 'Interpreter', 'latex');
ylabel('$$V(x)$$','Interpreter', 'latex');
grid on;
hold off;
- Solving and Plotting the Duffing Equation
% Define the system of ODEs for the non-chaotic Duffing equation
duffing_ode = @(t, X) [X(2);
X(1) - X(1).^3];
% Time span for the simulation
tspan = [0 10];
% Initial conditions [x(0), v(0)]
initial_conditions = [1; 1];
% Solve the ODE using ode45
[t, X] = ode45(duffing_ode, tspan, initial_conditions);
% Extract displacement (x) and velocity (v)
x = X(:, 1);
v = X(:, 2);
% Plot both x(t) and v(t) in the same figure
figure;
plot(t, x, 'b-', 'LineWidth', 2); % Plot x(t) with blue line
hold on;
plot(t, v, 'r--', 'LineWidth', 2); % Plot v(t) with red dashed line
% Add title, labels, and legend
title(' Component curve solutions to $$\ddot{x}-x+x^3=0$$','Interpreter', 'latex');
xlabel('t','Interpreter', 'latex');
ylabel('$$x(t) $$ and $$v(t) $$','Interpreter', 'latex');
legend('$$x(t)$$', ' $$v(t)$$','Interpreter', 'latex');
grid on;
hold off;
% Phase portrait with nullclines, equilibria, and vector field
figure;
hold on;
% Plot phase portrait
plot(x, v,'r', 'LineWidth', 2);
% Plot equilibrium points
plot([0 1 -1], [0 0 0], 'ro', 'MarkerSize', 5, 'MarkerFaceColor', 'g');
% Create a grid of points for the vector field
[x_vals, v_vals] = meshgrid(linspace(-2, 2, 20), linspace(-1, 1, 20));
% Compute the vector field components
dxdt = v_vals;
dvdt = x_vals - x_vals.^3;
% Plot the vector field
quiver(x_vals, v_vals, dxdt, dvdt, 'b');
% Set axis limits to [-1, 1]
xlim([-1.7 1.7]);
ylim([-1 1]);
% Labels and title
title('Phase-Plane solutions to $$\ddot{x}-x+x^3=0$$','Interpreter', 'latex');
xlabel('$$ (x)$$','Interpreter', 'latex');
ylabel('$$v(v)$$','Interpreter', 'latex');
grid on;
hold off;
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.
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
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
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]
% 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]
% 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');
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")
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")
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))
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")
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))
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.
Chord diagrams are very common in Python and R, but there are no related functions in MATLAB before. It is not easy to draw chord diagrams of the same quality as R language, But I created a MATLAB tool that could almost do it.
↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓
↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑
Here is the help document:
1 Data Format
The data requirement is a numerical matrix with all values greater than or equal to 0, or a table array, or a numerical matrix and cell array for names. First, give an example of a numerical matrix:
1.1 Numerical Matrix
dataMat=randi([0,5],[5,4]);
% 绘图(draw)
CC=chordChart(dataMat);
CC=CC.draw();
Since each object is not named, it will be automatically named Rn and Cn
1.2 Numerical Matrix and Cell Array for Names
dataMat=[2 0 1 2 5 1 2;
3 5 1 4 2 0 1;
4 0 5 5 2 4 3];
colName={'G1','G2','G3','G4','G5','G6','G7'};
rowName={'S1','S2','S3'};
CC=chordChart(dataMat,'rowName',rowName,'colName',colName);
CC=CC.draw();
RowName should be the same size as the rows of the matrix
ColName should be the same size as the columns of the matrix
For this example, if the value in the second row and third column is 1, it indicates that there is an energy flow from S2 to G3, and a chord with a width of 1 is needed between these two.
1.3 Table Array
A table array in the following format is required:
2 Decorate Chord
2.1 Batch modification of chords
Batch modification of chords can be done using the setChordProp function, and all properties of the Patch object can be modified. For example, modifying the color of the string, edge color, edge line sstyle, etc.:
CC.setChordProp('EdgeColor',[.3,.3,.3],'LineStyle','--',...
'LineWidth',.1,'FaceColor',[.3,.3,.3])
2.2 Individual Modification of Chord
The individual modification of chord can be done using the setChordMN function, where the values of m and n correspond exactly to the rows and columns of the original numerical matrix. For example, changing the color of the strings flowing from S2 to G4 to red:
CC.setChordMN(2,4,'FaceColor',[1,0,0])
2.3 Color Mapping of Chords
Just use function colormap to do so:
% version 1.7.0更新
% 可使用colormap函数直接修改颜色
% Colors can be adjusted directly using the function colormap(demo4)
colormap(flipud(pink))
3 Arc Shaped Block Decoration
3.1 Batch Decoration of Arc-Shaped Blocks
use:
- setSquareT_Prop
- setSquareF_Prop
to modify the upper and lower blocks separately, and all attributes of the Patch object can be modified. For example, batch modify the upper blocks (change to black):
CC.setSquareT_Prop('FaceColor',[0,0,0])
3.2 Arc-Shaped Blocks Individually Decoration
use:
- setSquareT_N
- setSquareF_N
to modify the upper and lower blocks separately. For example, modify the second block above separately (changed to red):
CC.setSquareT_N(2,'FaceColor',[.8,0,0])
4 Font Adjustment
Use the setFont function to adjust the font, and all properties of the text object can be modified. For example, changing the font size, font, and color of the text:
CC.setFont('FontSize',25,'FontName','Cambria','Color',[0,0,.8])
5 Show and Hide Ticks
Usage:
CC.tickState('on')
% CC.tickState('off')
6 Attribute 'Sep' with Adjustable Square Spacing
If the matrix size is large, the drawing will be out of scale:
dataMat=randi([0,1],[20,10]);
CC=chordChart(dataMat);
CC=CC.draw();
% CC.tickState('on')
We can modify its Sep attribute:
dataMat=randi([0,1],[20,10]);
% use Sep to decrease space (separation)
% 使用 sep 减小空隙
CC=chordChart(dataMat,'Sep',1/120);
CC=CC.draw();
7 Modify Text Direction
dataMat=randi([0,1],[20,10]);
% use Sep to decrease space (separation)
% 使用 sep 减小空隙
CC=chordChart(dataMat,'Sep',1/120);
CC=CC.draw();
CC.tickState('on')
% version 1.7.0更新
% 函数labelRatato用来旋转标签
% The function labelRatato is used to rotate the label
CC.labelRotate('on')
8 Add Tick Labels
dataMat=[2 0 1 2 5 1 2;
3 5 1 4 2 0 1;
4 0 5 5 2 4 3];
colName={'G1','G2','G3','G4','G5','G6','G7'};
rowName={'S1','S2','S3'};
CC=chordChart(dataMat,'rowName',rowName,'colName',colName);
CC=CC.draw();
CC.setFont('FontSize',17,'FontName','Cambria')
% 显示刻度和数值
% Displays scales and numeric values
CC.tickState('on')
CC.tickLabelState('on')
% 调节标签半径
% Adjustable Label radius
CC.setLabelRadius(1.3);
% figure()
% dataMat=[2 0 1 2 5 1 2;
% 3 5 1 4 2 0 1;
% 4 0 5 5 2 4 3];
% dataMat=dataMat+rand(3,7);
% dataMat(dataMat<1)=0;
%
% CC=chordChart(dataMat,'rowName',rowName,'colName',colName);
% CC=CC.draw();
% CC.setFont('FontSize',17,'FontName','Cambria')
%
% % 显示刻度和数值
% % Displays scales and numeric values
% CC.tickState('on')
% CC.tickLabelState('on')
%
% % 调节标签半径
% % Adjustable Label radius
% CC.setLabelRadius(1.4);
9 Custom Tick Label Format
A function handle is required to input numeric output strings. The format can be set through the setTickLabelFormat function, such as Scientific notation:
dataMat=[2 0 1 2 5 1 2;
3 5 1 4 2 0 1;
4 0 5 5 2 4 3];
dataMat=dataMat+rand(3,7);
dataMat(dataMat<1)=0;
dataMat=dataMat.*1000;
CC=chordChart(dataMat);
CC=CC.draw();
CC.setFont('FontSize',17,'FontName','Cambria')
% 显示刻度和数值
% Displays scales and numeric values
CC.tickState('on')
CC.tickLabelState('on')
% 调节标签半径
% Adjustable Label radius
CC.setLabelRadius(1.4);
% 调整数值字符串格式
% Adjust numeric string format
CC.setTickLabelFormat(@(x)sprintf('%0.1e',x))
10 A Demo
rng(2)
dataMat=randi([1,7],[11,5]);
colName={'Fly','Beetle','Leaf','Soil','Waxberry'};
rowName={'Bartomella','Bradyrhizobium','Dysgomonas','Enterococcus',...
'Lactococcus','norank','others','Pseudomonas','uncultured',...
'Vibrionimonas','Wolbachia'};
CC=chordChart(dataMat,'rowName',rowName,'colName',colName,'Sep',1/80);
CC=CC.draw();
% 修改上方方块颜色(Modify the color of the blocks above)
CListT=[0.7765 0.8118 0.5216;0.4431 0.4706 0.3843;0.5804 0.2275 0.4549;
0.4471 0.4039 0.6745;0.0157 0 0 ];
for i=1:5
CC.setSquareT_N(i,'FaceColor',CListT(i,:))
end
% 修改下方方块颜色(Modify the color of the blocks below)
CListF=[0.5843 0.6863 0.7843;0.1098 0.1647 0.3255;0.0902 0.1608 0.5373;
0.6314 0.7961 0.2118;0.0392 0.2078 0.1059;0.0157 0 0 ;
0.8549 0.9294 0.8745;0.3882 0.3255 0.4078;0.5020 0.7216 0.3843;
0.0902 0.1843 0.1804;0.8196 0.2314 0.0706];
for i=1:11
CC.setSquareF_N(i,'FaceColor',CListF(i,:))
end
% 修改弦颜色(Modify chord color)
for i=1:5
for j=1:11
CC.setChordMN(j,i,'FaceColor',CListT(i,:),'FaceAlpha',.5)
end
end
CC.tickState('on')
CC.labelRotate('on')
CC.setFont('FontSize',17,'FontName','Cambria')
Hope to have your Reviews and Stars!!!
↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓
↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑
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
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
% 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
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
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');
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, []);
How to leave feedback on a doc page
Leaving feedback is a two-step process. At the bottom of most pages in the MATLAB documentation is a star rating.
Start by selecting a star that best answers the question. After selecting a star rating, an edit box appears where you can offer specific feedback.
When you press "Submit" you'll see the confirmation dialog below. You cannot go back and edit your content, although you can refresh the page to go through that process again.
Tips on leaving feedback
- Be productive. The reader should clearly understand what action you'd like to see, what was unclear, what you think needs work, or what areas were really helpful.
- Positive feedback is also helpful. By nature, feedback often focuses on suggestions for changes but it also helps to know what was clear and what worked well.
- Point to specific areas of the page. This helps the reader to narrow the focus of the page to the area described by your feedback.
What happens to that feedback?
Before working at MathWorks I often left feedback on documentation pages but I never knew what happens after that. One day in 2021 I shared my speculation on the process:
> That feedback is received by MathWorks Gnomes which are never seen nor heard but visit the MathWorks documentation team at night while they are sleeping and whisper selected suggestions into their ears to manipulate their dreams. Occassionally this causes them to wake up with a Eureka moment that leads to changes in the documentation.
I'd like to let you in on the secret which is much less fanciful. Feedback left in the star rating and edit box are collected and periodically reviewed by the doc writers who look for trends on highly trafficked pages and finer grain feedback on less visited pages. Your feedback is important and often results in improvements.
A colleague said that you can search the Help Center using the phrase 'Introduced in' followed by a release version. Such as, 'Introduced in R2022a'. Doing this yeilds search results specific for that release.
Seems pretty handy so I thought I'd share.
Bringing the beauty of MathWorks Natick's tulips to life through code!
Remix challenge: create and share with us your new breeds of MATLAB tulips!
This cheat sheet is here:
reference:
- https://github.com/peijin94/matlabPlotCheatsheet
- https://github.com/mathworks/visualization-cheat-sheet
- https://www.mathworks.com/products/matlab/plot-gallery.html
- https://www.mathworks.com/help/matlab/release-notes.html
MATLAB used to have official visualization-cheat-sheet, but there have been quite a few new updates in MATLAB versions recently. Therefore, I made my own cheat sheet and marked the versions of each new thing that were released :
From Alpha Vantage's website: API Documentation | Alpha Vantage
Try using the built-in Matlab function webread(URL)... for example:
% copy a URL from the examples on the site
URL = 'https://www.alphavantage.co/query?function=TIME_SERIES_DAILY&symbol=IBM&apikey=demo'
% or use the pattern to create one
tickers = [{'IBM'} {'SPY'} {'DJI'} {'QQQ'}]; i = 1;
URL = ...
['https://www.alphavantage.co/query?function=TIME_SERIES_DAILY_ADJUSTED&outputsize=full&symbol=', ...
+ tickers{i}, ...
+ '&apikey=***Put Your API Key here***'];
X = webread(URL);
You can access any of the data available on the site as per the Alpha Vantage documentation using these two lines of code but with different designations for the requested data as per the documentation.
It's fun!
How to create a legend as follows?
Principle Explanation - Graphic Objects
Hidden Properties of Legend are laid as follows
In most cases, legends are drawn using LineLoop and Quadrilateral:
Both of these basic graphic objects are drawn in groups of four points, and the general principle is as follows:
Of course, you can arrange the points in order, or set VertexIndices whitch means the vertex order to obtain the desired quadrilateral shape:
Other objects
Compared to objects that can only be grouped into four points, we also need to introduce more flexible objects. Firstly, LineStrip, a graphical object that draws lines in groups of two points:
And TriangleStrip is a set of three points that draw objects to fill triangles, for example, complex polygons can be filled with multiple triangles:
Principle Explanation - Create and Replace
Let's talk about how to construct basic graphic objects, which are all constructed using undisclosed and very low-level functions, such as LineStrip, not through:
- LineStrip()
It is built through:
- matlab.graphics.primitive.world.LineStrip()
After building the object, the following properties must be set to make the hidden object visible:
- Layer
- ColorBinding
- ColorData
- VertexData
- PickableParts
The settings of these properties can refer to the original legend to form the object, which will not be elaborated here. You can also refer to the code I wrote.
Afterwards, set the newly constructed object's parent class as the Group parent class of the original component, and then hide the original component
newBoxEdgeHdl.Parent=oriBoxEdgeHdl.Parent;
oriBoxEdgeHdl.Visible='off';
The above is the entire process of component replacement, with two example codes written:
Semi transparent legend
function SPrettyLegend(lgd)
% Semitransparent rounded rectangle legend
% Copyright (c) 2023, Zhaoxu Liu / slandarer
% -------------------------------------------------------------------------
% Zhaoxu Liu / slandarer (2023). pretty legend
% (https://www.mathworks.com/matlabcentral/fileexchange/132128-pretty-legend),
% MATLAB Central File Exchange. 检索来源 2023/7/9.
% =========================================================================
if nargin<1
ax = gca;
lgd = get(ax,'Legend');
end
pause(1e-6)
Ratio = .1;
t1 = linspace(0,pi/2,4); t1 = t1([1,2,2,3,3,4]);
t2 = linspace(pi/2,pi,4); t2 = t2([1,2,2,3,3,4]);
t3 = linspace(pi,3*pi/2,4); t3 = t3([1,2,2,3,3,4]);
t4 = linspace(3*pi/2,2*pi,4); t4 = t4([1,2,2,3,3,4]);
XX = [1,1,1-Ratio+cos(t1).*Ratio,1-Ratio,Ratio,Ratio+cos(t2).*Ratio,...
0,0,Ratio+cos(t3).*Ratio,Ratio,1-Ratio,1-Ratio+cos(t4).*Ratio];
YY = [Ratio,1-Ratio,1-Ratio+sin(t1).*Ratio,1,1,1-Ratio+sin(t2).*Ratio,...
1-Ratio,Ratio,Ratio+sin(t3).*Ratio,0,0,Ratio+sin(t4).*Ratio];
% 圆角边框(border-radius)
oriBoxEdgeHdl = lgd.BoxEdge;
newBoxEdgeHdl = matlab.graphics.primitive.world.LineStrip();
newBoxEdgeHdl.AlignVertexCenters = 'off';
newBoxEdgeHdl.Layer = 'front';
newBoxEdgeHdl.ColorBinding = 'object';
newBoxEdgeHdl.LineWidth = 1;
newBoxEdgeHdl.LineJoin = 'miter';
newBoxEdgeHdl.WideLineRenderingHint = 'software';
newBoxEdgeHdl.ColorData = uint8([38;38;38;0]);
newBoxEdgeHdl.VertexData = single([XX;YY;XX.*0]);
newBoxEdgeHdl.Parent=oriBoxEdgeHdl.Parent;
oriBoxEdgeHdl.Visible='off';
% 半透明圆角背景(Semitransparent rounded background)
oriBoxFaceHdl = lgd.BoxFace;
newBoxFaceHdl = matlab.graphics.primitive.world.TriangleStrip();
Ind = [1:(length(XX)-1);ones(1,length(XX)-1).*(length(XX)+1);2:length(XX)];
Ind = Ind(:).';
newBoxFaceHdl.PickableParts = 'all';
newBoxFaceHdl.Layer = 'back';
newBoxFaceHdl.ColorBinding = 'object';
newBoxFaceHdl.ColorType = 'truecoloralpha';
newBoxFaceHdl.ColorData = uint8(255*[1;1;1;.6]);
newBoxFaceHdl.VertexData = single([XX,.5;YY,.5;XX.*0,0]);
newBoxFaceHdl.VertexIndices = uint32(Ind);
newBoxFaceHdl.Parent = oriBoxFaceHdl.Parent;
oriBoxFaceHdl.Visible = 'off';
end
Usage examples
clc; clear; close all
rng(12)
% 生成随机点(Generate random points)
mu = [2 3; 6 7; 8 9];
S = cat(3,[1 0; 0 2],[1 0; 0 2],[1 0; 0 1]);
r1 = abs(mvnrnd(mu(1,:),S(:,:,1),100));
r2 = abs(mvnrnd(mu(2,:),S(:,:,2),100));
r3 = abs(mvnrnd(mu(3,:),S(:,:,3),100));
% 绘制散点图(Draw scatter chart)
hold on
propCell = {'LineWidth',1.2,'MarkerEdgeColor',[.3,.3,.3],'SizeData',60};
scatter(r1(:,1),r1(:,2),'filled','CData',[0.40 0.76 0.60],propCell{:});
scatter(r2(:,1),r2(:,2),'filled','CData',[0.99 0.55 0.38],propCell{:});
scatter(r3(:,1),r3(:,2),'filled','CData',[0.55 0.63 0.80],propCell{:});
% 增添图例(Draw legend)
lgd = legend('scatter1','scatter2','scatter3');
lgd.Location = 'northwest';
lgd.FontSize = 14;
% 坐标区域基础修饰(Axes basic decoration)
ax=gca; grid on
ax.FontName = 'Cambria';
ax.Color = [0.9,0.9,0.9];
ax.Box = 'off';
ax.TickDir = 'out';
ax.GridColor = [1 1 1];
ax.GridAlpha = 1;
ax.LineWidth = 1;
ax.XColor = [0.2,0.2,0.2];
ax.YColor = [0.2,0.2,0.2];
ax.TickLength = [0.015 0.025];
% 隐藏轴线(Hide XY-Ruler)
pause(1e-6)
ax.XRuler.Axle.LineStyle = 'none';
ax.YRuler.Axle.LineStyle = 'none';
SPrettyLegend(lgd)
Heart shaped legend (exclusive to pie charts)
function pie2HeartLegend(lgd)
% Heart shaped legend for pie chart
% Copyright (c) 2023, Zhaoxu Liu / slandarer
% -------------------------------------------------------------------------
% Zhaoxu Liu / slandarer (2023). pretty legend
% (https://www.mathworks.com/matlabcentral/fileexchange/132128-pretty-legend),
% MATLAB Central File Exchange. 检索来源 2023/7/9.
% =========================================================================
if nargin<1
ax = gca;
lgd = get(ax,'Legend');
end
pause(1e-6)
% 心形曲线(Heart curve)
x = -1:1/100:1;
y1 = 0.6 * abs(x) .^ 0.5 + ((1 - x .^ 2) / 2) .^ 0.5;
y2 = 0.6 * abs(x) .^ 0.5 - ((1 - x .^ 2) / 2) .^ 0.5;
XX = [x, flip(x),x(1)]./3.4+.5;
YY = ([y1, y2,y1(1)]-.2)./2+.5;
Ind = [1:(length(XX)-1);2:length(XX)];
Ind = Ind(:).';
% 获取图例图标(Get Legend Icon)
lgdEntryChild = lgd.EntryContainer.NodeChildren;
iconSet = arrayfun(@(lgdEntryChild)lgdEntryChild.Icon.Transform.Children.Children,lgdEntryChild,UniformOutput=false);
% 基础边框句柄(Base Border Handle)
newEdgeHdl = matlab.graphics.primitive.world.LineStrip();
newEdgeHdl.AlignVertexCenters = 'off';
newEdgeHdl.Layer = 'front';
newEdgeHdl.ColorBinding = 'object';
newEdgeHdl.LineWidth = .8;
newEdgeHdl.LineJoin = 'miter';
newEdgeHdl.WideLineRenderingHint = 'software';
newEdgeHdl.ColorData = uint8([38;38;38;0]);
newEdgeHdl.VertexData = single([XX;YY;XX.*0]);
newEdgeHdl.VertexIndices = uint32(Ind);
% 基础多边形面句柄(Base Patch Handle)
newFaceHdl = matlab.graphics.primitive.world.TriangleStrip();
Ind = [1:(length(XX)-1);ones(1,length(XX)-1).*(length(XX)+1);2:length(XX)];
Ind = Ind(:).';
newFaceHdl.PickableParts = 'all';
newFaceHdl.Layer = 'middle';
newFaceHdl.ColorBinding = 'object';
newFaceHdl.ColorType = 'truecoloralpha';
newFaceHdl.ColorData = uint8(255*[1;1;1;.6]);
newFaceHdl.VertexData = single([XX,.5;YY,.5;XX.*0,0]);
newFaceHdl.VertexIndices = uint32(Ind);
% 替换图例图标(Replace Legend Icon)
for i = 1:length(iconSet)
oriEdgeHdl = iconSet{i}(1);
tNewEdgeHdl = copy(newEdgeHdl);
tNewEdgeHdl.ColorData = oriEdgeHdl.ColorData;
tNewEdgeHdl.Parent = oriEdgeHdl.Parent;
oriEdgeHdl.Visible = 'off';
oriFaceHdl = iconSet{i}(2);
tNewFaceHdl = copy(newFaceHdl);
tNewFaceHdl.ColorData = oriFaceHdl.ColorData;
tNewFaceHdl.Parent = oriFaceHdl.Parent;
oriFaceHdl.Visible = 'off';
end
end
Usage examples
clc; clear; close all
% 生成随机点(Generate random points)
X = [1 3 0.5 2.5 2];
pieHdl = pie(X);
% 修饰饼状图(Decorate pie chart)
colorList=[0.4941 0.5490 0.4118
0.9059 0.6510 0.3333
0.8980 0.6157 0.4980
0.8902 0.5137 0.4667
0.4275 0.2824 0.2784];
for i = 1:2:length(pieHdl)
pieHdl(i).FaceColor=colorList((i+1)/2,:);
pieHdl(i).EdgeColor=colorList((i+1)/2,:);
pieHdl(i).LineWidth=1;
pieHdl(i).FaceAlpha=.6;
end
for i = 2:2:length(pieHdl)
pieHdl(i).FontSize=13;
pieHdl(i).FontName='Times New Roman';
end
lgd=legend('FontSize',13,'FontName','Times New Roman','TextColor',[1,1,1].*.3);
pie2HeartLegend(lgd)
Temporary print statements are often helpful during debugging but it's easy to forget to remove the statements or sometimes you may not have writing privileges for the file. This tip uses conditional breakpoints to add print statements without ever editing the file!
What are conditional breakpoints?
Conditional breakpoints allow you to write a conditional statement that is executed when the selected line is hit and if the condition returns true, MATLAB pauses at that line. Otherwise, it continues.
The Hack: use ~fprintf() as the condition
fprintf prints information to the command window and returns the size of the message in bytes. The message size will always be greater than 0 which will always evaluate as true when converted to logical. Therefore, by negating an fprintf statement within a conditional breakpoint, the fprintf command will execute, print to the command window, and evalute as false which means the execution will continue uninterupted!
How to set a conditional break point
1. Right click the line number where you want the condition to be evaluated and select "Set Conditional Breakpoint"
2. Enter a valid MATLAB expression that returns a logical scalar value in the editor dialog.
Handy one-liners
Check if a line is reached: Don't forget the negation (~) and the line break (\n)!
~fprintf('Entered callback function\n')
Display the call stack from the break point line: one of my favorites!
~fprintf('%s\n',formattedDisplayText(struct2table(dbstack)))
Inspect variable values: For scalar values,
~fprintf('v = %.5f\n', v)
~fprintf('%s\n', formattedDisplayText(v)).
Make sense of frequent hits: In some situations such as responses to listeners or interactive callbacks, a line can be executed 100s of times per second. Incorporate a timestamp to differentiate messages during rapid execution.
~fprintf('WindowButtonDownFcn - %s\n', datetime('now'))
Closing
This tip not only keeps your code clean but also offers a dynamic way to monitor code execution and variable states without permanent modifications. Interested in digging deeper? @Steve Eddins takes this tip to the next level with his Code Trace for MATLAB tool available on the File Exchange (read more).
Summary animation
To reproduce the events in this animation:
% buttonDownFcnDemo.m
fig = figure();
tcl = tiledlayout(4,4,'TileSpacing','compact');
for i = 1:16
ax = nexttile(tcl);
title(ax,"#"+string(i))
ax.ButtonDownFcn = @axesButtonDownFcn;
xlim(ax,[-1 1])
ylim(ax,[-1,1])
hold(ax,'on')
end
function axesButtonDownFcn(obj,event)
colors = lines(16);
plot(obj,event.IntersectionPoint(1),event.IntersectionPoint(2),...
'ko','MarkerFaceColor',colors(obj.Layout.Tile,:))
end
If you've dabbled in "procedural generation," (algorithmically generating natural features), you may have come across the problem of sphere texturing. How to seamlessly texture a sphere is not immediately obvious. Watch what happens, for example, if you try adding power law noise to an evenly sampled grid of spherical angle coordinates (i.e. a "UV sphere" in Blender-speak):
% Example: how [not] to texture a sphere:
rng(2, 'twister'); % Make what I have here repeatable for you
% Make our radial noise, mapped onto an equal spaced longitude and latitude
% grid.
N = 51;
b = linspace(-1, 1, N).^2;
r = abs(ifft2(exp(6i*rand(N))./(b'+b+1e-5))); % Power law noise
r = rescale(r, 0, 1) + 5;
[lon, lat] = meshgrid(linspace(0, 2*pi, N), linspace(-pi/2, pi/2, N));
[x2, y2, z2] = sph2cart(lon, lat, r);
r2d = @(x)x*180/pi;
% Radial surface texture
subplot(1, 3, 1);
imagesc(r, 'Xdata', r2d(lon(1,:)), 'Ydata', r2d(lat(:, 1)));
xlabel('Longitude (Deg)');
ylabel('Latitude (Deg)');
title('Texture (radial variation)');
% View from z axis
subplot(1, 3, 2);
surf(x2, y2, z2, r);
axis equal
view([0, 90]);
title('Top view');
% Side view
subplot(1, 3, 3);
surf(x2, y2, z2, r);
axis equal
view([-90, 0]);
title('Side view');
The created surface shows "pinching" at the poles due to different radial values mapping to the same location. Furthermore, the noise statistics change based on the density of the sampling on the surface.
How can this be avoided? One standard method is to create a textured volume and sample the volume at points on a sphere. Code for doing this is quite simple:
rng default % Make our noise realization repeatable
% Create our 3D power-law noise
N = 201;
b = linspace(-1, 1, N);
[x3, y3, z3] = meshgrid(b, b, b);
b3 = x3.^2 + y3.^2 + z3.^2;
r = abs(ifftn(ifftshift(exp(6i*randn(size(b3)))./(b3.^1.2 + 1e-6))));
% Modify it - make it more interesting
r = rescale(r);
r = r./(abs(r - 0.5) + .1);
% Sample on a sphere
[x, y, z] = sphere(500);
% Plot
ir = interp3(x3, y3, z3, r, x, y, z, 'linear', 0);
surf(x, y, z, ir);
shading flat
axis equal off
set(gcf, 'color', 'k');
colormap(gray);
The result of evaluating this code is a seamless, textured sphere with no discontinuities at the poles or variation in the spatial statistics of the noise texture:
But what if you want to smooth it or perform some other local texture modification? Smoothing the volume and resampling is not equivalent to smoothing the surficial features shown on the map above.
A more flexible alternative is to treat the samples on the sphere surface as a set of interconnected nodes that are influenced by adjacent values. Using this approach we can start by defining the set of nodes on a sphere surface. These can be sampled almost arbitrarily, though the noise statistics will vary depending on the sampling strategy.
One noise realisation I find attractive can be had by randomly sampling a sphere. Normalizing a point in N-dimensional space by its 2-norm projects it to the surface of an N-dimensional unit sphere, so randomly sampling a sphere can be done very easily using randn() and vecnorm():
N = 5e3; % Number of nodes on our sphere
g=randn(3,N); % Random 3D points around origin
p=g./vecnorm(g); % Projected to unit sphere
The next step is to find each point's "neighbors." The first step is to find the convex hull. Since each point is on the sphere, the convex hull will include each point as a vertex in the triangulation:
k=convhull(p');
In the above, k is an N x 3 set of indices where each row represents a unique triangle formed by a triplicate of points on the sphere surface. The vertices of the full set of triangles containing a point describe the list of neighbors to that point.
What we want now is a large, sparse symmetric matrix where the indices of the columns & rows represent the indices of the points on the sphere and the nth row (and/or column) contains non-zero entries at the indices corresponding to the neighbors of the nth point.
How to do this? You could set up a tiresome nested for-loop searching for all rows (triangles) in k that contain some index n, or you could directly index via:
c=@(x)sparse(k(:,x)*[1,1,1],k,1,N,N);
t=c(1)|c(2)|c(3);
The result is the desired sparse connectivity matrix: a matrix with non-zero entries defining neighboring points.
So how do we create a textured sphere with this connectivity matrix? We will use it to form a set of equations that, when combined with the concept of "regularization," will allow us to determine the properties of the randomness on the surface. Our regularizer will penalize the difference of the radial distance of a point and the average of its neighbors. To do this we replace the main diagonal with the negative of the sum of the off-diagonal components so that the rows and columns are zero-mean. This can be done via:
w=spdiags(-sum(t,2)+1,0,double(t));
Now we invoke a bit of linear algebra. Pretend x is an N-length vector representing the radial distance of each point on our sphere with the noise realization we desire. Y will be an N-length vector of "observations" we are going to generate randomly, in this case using a uniform distribution (because it has a bias and we want a non-zero average radius, but you can play around with different distributions than uniform to get different effects):
Y=rand(N,1);
and A is going to be our "transformation" matrix mapping x to our noisy observations:
Ax = Y
In this case both x and Y are N length vectors and A is just the identity matrix:
A = speye(N);
Y, however, doesn't create the noise realization we want. So in the equation above, when solving for x we are going to introduce a regularizer which is going to penalize unwanted behavior of x by some amount. That behavior is defined by the point-neighbor radial differences represented in matrix w. Our estimate of x can then be found using one of my favorite Matlab assets, the "\" operator:
smoothness = 10; % Smoothness penalty: higher is smoother
x = (A+smoothness*w'*w)\Y; % Solving for radii
The vector x now contains the radii with the specified noise realization for the sphere which can be created simply by multiplying x by p and plotting using trisurf:
p2 = p.*x';
trisurf(k,p2(1,:),p2(2,:),p2(3,:),'FaceC', 'w', 'EdgeC', 'none','AmbientS',0,'DiffuseS',0.6,'SpecularS',1);
light;
set(gca, 'color', 'k');
axis equal
The following images show what happens as you change the smoothness parameter using values [.1, 1, 10, 100] (left to right):
Now you know a couple ways to make a textured sphere: that's the starting point for having a lot of fun with basic procedural planet, moon, or astroid generation! Here's some examples of things you can create based on these general ideas:
I am often talking to new MATLAB users. I have put together one script. If you know how this script works, why, and what each line means, you will be well on your way on your MATLAB learning journey.
% Clear existing variables and close figures
clear;
close all;
% Print to the Command Window
disp('Hello, welcome to MATLAB!');
% Create a simple vector and matrix
vector = [1, 2, 3, 4, 5];
matrix = [1, 2, 3; 4, 5, 6; 7, 8, 9];
% Display the created vector and matrix
disp('Created vector:');
disp(vector);
disp('Created matrix:');
disp(matrix);
% Perform element-wise multiplication
result = vector .* 2;
% Display the result of the operation
disp('Result of element-wise multiplication of the vector by 2:');
disp(result);
% Create plot
x = 0:0.1:2*pi; % Generate values from 0 to 2*pi
y = sin(x); % Calculate the sine of these values
% Plotting
figure; % Create a new figure window
plot(x, y); % Plot x vs. y
title('Simple Plot of sin(x)'); % Give the plot a title
xlabel('x'); % Label the x-axis
ylabel('sin(x)'); % Label the y-axis
grid on; % Turn on the grid
disp('This is the end of the script. Explore MATLAB further to learn more!');
Introduction
Comma-separated lists are really very simple. You use them all the time. Here is one:
a,b,c,d
That is a comma-separated list containing four variables, the variables a, b, c, and d. Every time you write a list separated by commas then you are writing a comma-separated list. Most commonly you would write a comma-separated list as inputs when calling a function:
fun(a,b,c,d)
or as arguments to the concatenation operator or cell construction operator:
[a,b,c,d]
{a,b,c,d}
or as function outputs:
[a,b,c,d] = fun();
It is very important to understand that in general a comma-separated list is NOT one variable (but it could be). However, sometimes it is useful to create a comma-separated list from one variable (or define one variable from a comma-separated list), and MATLAB has several ways of doing this from various container array types:
struct_array.field % all elements
struct_array(idx).field % selected elements
cell_array{:} % all elements
cell_array{idx} % selected elements
string_array{:} % all elements
string_array{idx} % selected elements
Note that in all cases, the comma-separated list consists of the content of the container array, not subsets (or "slices") of the container array itself (use parentheses to "slice" any array). In other words, they will be equivalent to writing this comma-separated list of the container array content:
content1, content2, content3, .. , contentN
and will return as many content arrays as the original container array has elements (or that you select using indexing, in the requested order). A comma-separated list of one element is just one array, but in general there can be any number of separate arrays in the comma-separated list (zero, one, two, three, four, or more). Here is an example showing that a comma-separated list generated from the content of a cell array is the same as a comma-separated list written explicitly:
>> C = {1,0,Inf};
>> C{:}
ans =
1
ans =
0
ans =
Inf
>> 1,0,Inf
ans =
1
ans =
0
ans =
Inf
How to Use Comma-Separated Lists
Function Inputs: Remember that every time you call a function with multiple input arguments you are using a comma-separated list:
fun(a,b,c,d)
and this is exactly why they are useful: because you can specify the arguments for a function or operator without knowing anything about the arguments (even how many there are). Using the example cell array from above:
>> vertcat(C{:})
ans =
1
0
Inf
which, as we should know by now, is exactly equivalent to writing the same comma-separated list directly into the function call:
>> vertcat(1,0,Inf)
ans =
1
0
Inf
How can we use this? Commonly these are used to generate vectors of values from a structure or cell array, e.g. to concatenate the filenames which are in the output structure of dir:
S = dir(..);
F = {S.name}
which is simply equivalent to
F = {S(1).name, S(2).name, S(3).name, .. , S(end).name}
Or, consider a function with multiple optional input arguments:
opt = {'HeaderLines',2, 'Delimiter',',', 'CollectOutputs',true);
fid = fopen(..);
C = textscan(fid,'%f%f',opt{:});
fclose(fid);
Note how we can pass the optional arguments as a comma-separated list. Remember how a comma-separated list is equivalent to writing var1,var2,var3,..., then the above example is really just this:
C = textscan(fid,'%f%f', 'HeaderLines',2, 'Delimiter',',', 'CollectOutputs',true)
with the added advantage that we can specify all of the optional arguments elsewhere and handle them as one cell array (e.g. as a function input, or at the top of the file). Or we could select which options we want simply by using indexing on that cell array. Note that varargin and varargout can also be useful here.
Function Outputs: In much the same way that the input arguments can be specified, so can an arbitrary number of output arguments. This is commonly used for functions which return a variable number of output arguments, specifically ind2sub and gradient and ndgrid. For example we can easily get all outputs of ndgrid, for any number of inputs (in this example three inputs and three outputs, determined by the number of elements in the cell array):
C = {1:3,4:7,8:9};
[C{:}] = ndgrid(C{:});
which is thus equivalent to:
[C{1},C{2},C{3}] = ndgrid(C{1},C{2},C{3});
Further Topics:
MATLAB documentation:
Click on these links to jump to relevant comments below:
Dynamic Indexing (indexing into arrays with arbitrary numbers of dimensions)
Summary
Just remember that in general a comma-separated list is not one variable (although they can be), and that they are exactly what they say: a list (of arrays) separated with commas. You use them all the time without even realizing it, every time you write this:
fun(a,b,c,d)
Puoi anche selezionare un sito web dal seguente elenco:
Americhe
- América Latina (Español)
- Canada (English)
- United States (English)
Europa
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia-Pacifico
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)