ディスカッションの検索

J.K043006
最後のアクティビティ: 2024 年 8 月 29 日

"corr" malfunctioning in cody

I've been working on some matrix problems recently(Problem 55225)

and this is my code

It turns out that "Undefined function 'corr' for input arguments of type 'double'." However, should't the input argument of "corr" be column vectors with single/double values? What's even going on there?

Chris Hooper
最後のアクティビティ: 2024 年 8 月 29 日

uniqueN function would be useful

isequaln exists to return true when NaN==NaN.

unique treats NaN==NaN as false (as it should) requiring NaN to be replaced if NaN is not considered unique in a particular application. In my application, I am checking uniqueness of table rows using [table_unique,index_unique]=unique(table,"rows","sorted") and would prefer to keep NaN as NaN or missing in table_unique without the overhead of replacing it with a dummy value then replacing it again. Dummy values also have the risk of matching existing values in the table, requiring first finding a dummy value that is not in the table.

uniquen (similar to isequaln) would be more eloquent.

Please point out if I am missing something!

Athanasios Paraskevopoulos
最後のアクティビティ: 2024 年 8 月 17 日

Duffing equation : Transition to Chaos

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

Then 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:

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;

Athanasios Paraskevopoulos
最後のアクティビティ: 2024 年 8 月 15 日

The Non-Chaotic Duffing Equation

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;

Athanasios Paraskevopoulos
最後のアクティビティ: 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

Off Topic

David
最後のアクティビティ: 2024 年 8 月 8 日

Visual PDE

A library of runnable PDEs. See the equations! Modify the parameters! Visualize the resulting system in your browser! Convenient, fast, and instructive.

Visit https://visualpde.com/explore.html

Athanasios Paraskevopoulos
最後のアクティビティ: 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
最後のアクティビティ: 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
最後のアクティビティ: 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.)

Off Topic

Athanasios Paraskevopoulos
最後のアクティビティ: 2024 年 12 月 10 日 6:40

Looking for Beginner Courses on Drone Programming

Hi everyone,

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

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

Thank you!

Off Topic

Muhammad
最後のアクティビティ: 2024 年 8 月 26 日

What do you like most about MATLAB?

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

DH
最後のアクティビティ: 2024 年 8 月 8 日

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

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

Marisa
最後のアクティビティ: 2024 年 8 月 26 日

Can't click Roll the Dice!

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

Gregory Vernon
最後のアクティビティ: 2024 年 7 月 13 日

More easily enable CAD-like view mode

Something that had bothered me ever since I became an FEA analyst (2012) was the apparent inability of the "camera" in Matlab's 3D plot to function like the "cameras" in CAD/CAE packages.

For instance, load the ForearmLink.stl model that ships with the PDE Toolbox in Matlab and ParaView and try rotating the model.

clear

close all

gm = importGeometry( "ForearmLink.stl" );

pdegplot(gm)

To provide talking points, here's a YouTube video I recorded.

Things to observe:

Note that I cant seem to rotate continuously around the x-axis. It appears to only support rotations from [0, 360] as opposed to [-inf, inf]. So, for example, if I'm looking in the Y+ direction and rotate around X so that I'm looking at the Z- direction, and then want to look in the Y- direction, I can't simply keep rotating around the X axis... instead have to rotate 180 degrees around the Z axis and then around the X axis. I'm not aware of any data visualization applications (e.g., ParaView, VisIt, EnSight) or CAD/CAE tools with such an interaction.
Note that at the 50 second mark, I set a view in ParaView: looking in the [X-, Y-, Z-] direction with Y+ up. Try as I might in Matlab, I'm unable to achieve that same view perspective.

Today I discovered that if one turns on the Camera Toolbar from the View menubar, then clicks the Orbit Camera icon, then the No Principal Axis icon:

That then it acts in the manner I've long desired. Oh, and also, for the interested, it is programmatically available: https://www.mathworks.com/help/matlab/ref/cameratoolbar.html

I might humbly propose this mode either be made more discoverable, similar to the little interaction widgets that pop up in figures:

Or maybe use the middle-mouse button to temporarily use this mode (a mouse setting in, e.g., Abaqus/CAE).

Honzik
最後のアクティビティ: 2024 年 7 月 18 日

Matlab Custom Font

I've noticed is that the highly rated fonts for coding (e.g. Fira Code, Inconsolata, etc.) seem to overlook one issue that is key for coding in Matlab. While these fonts make 0 and O, as well as the 1 and l easily distinguishable, the brackets are not. Quite often the curly bracket looks similar to the curved bracket, which can lead to mistakes when coding or reviewing code.

So I was thinking: Could Mathworks put together a team to review good programming fonts, and come up with their own custom font designed specifically and optimized for Matlab syntax?

jmgoldba
最後のアクティビティ: 2024 年 6 月 21 日

Add Tools->Basic Fitting->...Weighted least squares

An option for 10th degree polynomials but no weighted linear least squares. Seriously? Jesse

Athanasios Paraskevopoulos
最後のアクティビティ: 2024 年 7 月 4 日

Stochastic simulation of a novel pathogen

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

Simulating the stochastic model 100 times for

Since γ is 0.4 per day,

per day

% Define the parameters

beta = 0.36;

gamma = 0.4;

n_sims = 100;

tf = 100; % Time frame changed to 100

% Calculate R0

R0 = beta / gamma

R0 = 0.9000

% Initial state values

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

% Define the propensities and state change matrix

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

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

% Define the Gillespie algorithm function

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

t = 0;

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

t_values = t;

state_values = state';

while t < tf

rates = a(state);

rate_sum = sum(rates);

if rate_sum == 0

break;

end

tau = -log(rand) / rate_sum;

t = t + tau;

r = rand * rate_sum;

cum_sum_rates = cumsum(rates);

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

state = state + nu(:, reaction_index);

% Update cumulative incidence if infection occurred

if reaction_index == 1

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

end

t_values = [t_values; t];

state_values = [state_values; state'];

end

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

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

% Define the propensities and state change matrix

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

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

% Set random seed for reproducibility

rng(11);

% Initialize plot

figure;

hold on;

for i = 1:n_sims

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

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

while length(t) == 1

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

end

if strcmp(plot_type, 'cumulative_incidence')

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

elseif strcmp(plot_type, 'prevalence')

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

end

xlabel('Time (days)');

if strcmp(plot_type, 'cumulative_incidence')

ylabel('Cumulative Incidence');

ylim([0 inf]);

elseif strcmp(plot_type, 'prevalence')

ylabel('Prevalence of Infection');

ylim([0 50]);

end

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

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

xlim([0 tf]);

grid on;

hold off;

end

% Simulate the model 100 times and plot cumulative incidence

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

% Simulate the model 100 times and plot prevalence

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

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

Ned Gulley
最後のアクティビティ: 2024 年 6 月 13 日

Make a video: "How I Solved It"

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

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

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

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

goc3
最後のアクティビティ: 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.