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

Hi everyone, I wrote several fancy functions that may help your coding experience, since they are in very early developing stage, I will be thankful if anyone could try them and give some feedbacks. Currently I have following:

- fstr: a Python f-string like expression
- printf: an easy to use fprintf function, accepts multiple arguments with seperator, end string control.

I will bring more functions or packages like logger when I am available.

The code is open sourced at GitHub with simple examples: https://github.com/bentrainer/MMGA

MATLAB Comprehensive commands list:

- clc - clears command window, workspace not affected
- clear - clears all variables from workspace, all variable values are lost
- diary - records into a file almost everything that appears in command window.
- exit - exits the MATLAB session
- who - prints all variables in the workspace
- whos - prints all variables in current workspace, with additional information.

Ch. 2 - Basics:

- Mathematical constants: pi, i, j, Inf, Nan, realmin, realmax
- Precedence rules: (), ^, negation, */, +-, left-to-right
- and, or, not, xor
- exp - exponential
- log - natural logarithm
- log10 - common logarithm (base 10)
- sqrt (square root)

- fprintf("final amount is %f units.", amt);
- can have: %f, %d, %i, %c, %s
- %f - fixed-point notation
- %e - scientific notation with lowercase e
- disp - outputs to a command window
- % - fieldWith.precision convChar

- MyArray = [startValue : IncrementingValue : terminatingValue]

Linspace

linspace(xStart, xStop, numPoints)

% xStart: Starting value

% xStop: Stop value

% numPoints: Number of linear-spaced points, including xStart and xStop

% Outputs a row array with the resulting values

Logspace

logspace(powerStart, powerStop, numPoints)

% powerStart: Starting value 10^powerStart

% powerStop: Stop value 10^powerStop

% numPoints: Number of logarithmic spaced points, including 10^powerStart and 10^powerStop

% Outputs a row array with the resulting values

- Transpose an array with []'

Element-Wise Operations

rowVecA = [1, 4, 5, 2];

rowVecB = [1, 3, 0, 4];

sumEx = rowVecA + rowVecB % Element-wise addition

diffEx = rowVecA - rowVecB % Element-wise subtraction

dotMul = rowVecA .* rowVecB % Element-wise multiplication

dotDiv = rowVecA ./ rowVecB % Element-wise division

dotPow = rowVecA .^ rowVecB % Element-wise exponentiation

- isinf(A) - check if the array elements are infinity
- isnan(A)

Rounding Functions

- ceil(x) - rounds each element of x to nearest integer >= to element
- floor(x) - rounds each element of x to nearest integer <= to element
- fix(x) - rounds each element of x to nearest integer towards 0
- round(x) - rounds each element of x to nearest integer. if round(x, N), rounds N digits to the right of the decimal point.
- rem(dividend, divisor) - produces a result that is either 0 or has the same sign as the dividen.
- mod(dividend, divisor) - produces a result that is either 0 or same result as divisor.
- Ex: 12/2, 12 is dividend, 2 is divisor
- sum(inputArray) - sums all entires in array

Complex Number Functions

- abs(z) - absolute value, is magnitude of complex number (phasor form r*exp(j0)
- angle(z) - phase angle, corresponds to 0 in r*exp(j0)
- complex(a,b) - creates complex number z = a + jb
- conj(z) - given complex conjugate a - jb
- real(z) - extracts real part from z
- imag(z) - extracts imaginary part from z
- unwrap(z) - removes the modulus 2pi from an array of phase angles.

Statistics Functions

- mean(xAr) - arithmetic mean calculated.
- std(xAr) - calculated standard deviation
- median(xAr) - calculated the median of a list of numbers
- mode(xAr) - calculates the mode, value that appears most often
- max(xAr)
- min(xAr)
- If using &&, this means that if first false, don't bother evaluating second

Random Number Functions

- rand(numRand, 1) - creates column array
- rand(1, numRand) - creates row array, both with numRand elements, between 0 and 1
- randi(maxRandVal, numRan, 1) - creates a column array, with numRand elements, between 1 and maxRandValue.
- randn(numRand, 1) - creates a column array with normally distributed values.
- Ex: 10 * rand(1, 3) + 6
- "10*rand(1, 3)" produces a row array with 3 random numbers between 0 and 10. Adding 6 to each element results in random values between 6 and 16.
- randi(20, 1, 5)
- Generates 5 (last argument) random integers between 1 (second argument) and 20 (first argument). The arguments 1 and 5 specify a 1 × 5 row array is returned.

Discrete Integer Mathematics

- primes(maxVal) - returns prime numbers less than or equal to maxVal
- isprime(inputNums) - returns a logical array, indicating whether each element is a prime number
- factor(intVal) - returns the prime factors of a number
- gcd(aVals, bVals) - largest integer that divides both a and b without a remainder
- lcm(aVals, bVals) - smallest positive integer that is divisible by both a and b
- factorial(intVals) - returns the factorial
- perms(intVals) - returns all the permutations of the elements int he array intVals in a 2D array pMat.
- randperm(maxVal)
- nchoosek(n, k)
- binopdf(x, n, p)

Concatenation

- cat, vertcat, horzcat
- Flattening an array, becomes vertical: sampleList = sampleArray ( : )

Dimensional Properties of Arrays

- nLargest = length(inArray) - number of elements along largest dimension
- nDim = ndims(inArray)
- nArrElement = numel(inArray) - nuber of array elements
- [nRow, nCol] = size(inArray) - returns the number of rows and columns on array. use (inArray, 1) if only row, (inArray, 2) if only column needed
- aZero = zeros(m, n) - creates an m by n array with all elements 0
- aOnes = ones(m, n) - creates an m by n array with all elements set to 1
- aEye = eye(m, n) - creates an m by n array with main diagonal ones
- aDiag = diag(vector) - returns square array, with diagonal the same, 0s elsewhere.
- outFlipLR = fliplr(A) - Flips array left to right.
- outFlipUD = flipud(A) - Flips array upside down.
- outRot90 = rot90(A) - Rotates array by 90 degrees counter clockwise around element at index (1,1).
- outTril = tril(A) - Returns the lower triangular part of an array.
- outTriU = triu(A) - Returns the upper triangular part of an array.
- arrayOut = repmat(subarrayIn, mRow, nCol), creates a large array by replicating a smaller array, with mRow x nCol tiling of copies of subarrayIn
- reshapeOut - reshape(arrayIn, numRow, numCol) - returns array with modifid dimensions. Product must equal to arrayIn of numRow and numCol.
- outLin = find(inputAr) - locates all nonzero elements of inputAr and returns linear indices of these elements in outLin.
- [sortOut, sortIndices] = sort(inArray) - sorts elements in ascending order, results result in sortOut. specify 'descend' if you want descending order. sortIndices hold the sorted indices of the array elements, which are row indices of the elements of sortOut in the original array
- [sortOut, sortIndices] = sortrows(inArray, colRef) - sorts array based on values in column colRef while keeping the rows together. Bases on first column by default.
- isequal(inArray1, inarray2, ..., inArrayN)
- isequaln(inArray1, inarray2, ..., inarrayn)

- arrayChar = ischar(inArray) - ischar tests if the input is a character array.
- arrayLogical = islogical(inArray) - islogical tests for logical array.
- arrayNumeric = isnumeric(inArray) - isnumeric tests for numeric array.
- arrayInteger = isinteger(inArray) - isinteger tests whether the input is integer type (Ex: uint8, uint16, ...)
- arrayFloat = isfloat(inArray) - isfloat tests for floating-point array.
- arrayReal= isreal(inArray) - isreal tests for real array.
- objIsa = isa(obj,ClassName) - isa determines whether input obj is an object of specified class ClassName.
- arrayScalar = isscalar(inArray) - isscalar tests for scalar type.
- arrayVector = isvector(inArray) - isvector tests for a vector (a 1D row or column array).
- arrayColumn = iscolumn(inArray) - iscolumn tests for column 1D arrays.
- arrayMatrix = ismatrix(inArray) - ismatrix returns true for a scalar or array up to 2D, but false for an array of more than 2 dimensions.
- arrayEmpty = isempty(inArray) - isempty tests whether inArray is empty.
- primeArray = isprime(inArray) - isprime returns a logical array primeArray, of the same size as inArray. The value at primeArray(index) is true when inArray(index) is a prime number. Otherwise, the values are false.
- finiteArray = isfinite(inArray) - isfinite returns a logical array finiteArray, of the same size as inArray. The value at finiteArray(index) is true when inArray(index) is finite. Otherwise, the values are false.
- infiniteArray = isinf(inArray) - isinf returns a logical array infiniteArray, of the same size as inArray. The value at infiniteArray(index) is true when inArray(index) is infinite. Otherwise, the values are false.
- nanArray = isnan(inArray) - isnan returns a logical array nanArray, of the same size as inArray. The value at nanArray(index) is true when inArray(index) is NaN. Otherwise, the values are false.
- allNonzero = all(inArray) - all identifies whether all array elements are non-zero (true). Instead of testing elements along the columns, all(inArray, 2) tests along the rows. all(inArray,1) is equivalent to all(inArray).
- anyNonzero = any(inArray) - any identifies whether any array elements are non-zero (true), and false otherwise. Instead of testing elements along the columns, any(inArray, 2) tests along the rows. any(inArray,1) is equivalent to any(inArray).
- logicArray = ismember(inArraySet,areElementsMember) - ismember returns a logical array logicArray the same size as inArraySet. The values at logicArray(i) are true where the elements of the first array inArraySet are found in the second array areElementsMember. Otherwise, the values are false. Similar values found by ismember can be extracted with inArraySet(logicArray).
- any(x) - Returns true if x is nonzero; otherwise, returns false.
- isnan(x) - Returns true if x is NaN (Not-a-Number); otherwise, returns false.
- isfinite(x) - Returns true if x is finite; otherwise, returns false. Ex: isfinite(Inf) is false, and isfinite(10) is true.
- isinf(x) - Returns true if x is +Inf or -Inf; otherwise, returns false.

Relational Operators

a & b, and(a, b)

a | b, or(a, b)

~a, not(a)

xor(a, b)

- fctName = @(arglist) expression - anonymous function
- nargin - keyword returns the number of input arguments passed to the function.

Looping

while condition

% code

end

for index = startVal:endVal

% code

end

- continue: Skips the rest of the current loop iteration and begins the next iteration.
- break: Exits a loop before it has finished all iterations.

switch expression

case value1

% code

case value2

% code

otherwise

% code

end

Comprehensive Overview (may repeat)

Built in functions/constants

abs(x) - absolute value

pi - 3.1415...

inf - ∞

eps - floating point accuracy 1e6 106

sum(x) - sums elements in x

cumsum(x) - Cumulative sum

prod - Product of array elements cumprod(x) cumulative product

diff - Difference of elements round/ceil/fix/floor Standard functions..

*Standard functions: sqrt, log, exp, max, min, Bessel *Factorial(x) is only precise for x < 21

Variable Generation

j:k - row vector

j:i:k - row vector incrementing by i

linspace(a,b,n) - n points linearly spaced and including a and b

NaN(a,b) - axb matrix of NaN values

ones(a,b) - axb matrix with all 1 values

zeros(a,b) - axb matrix with all 0 values

meshgrid(x,y) - 2d grid of x and y vectors

global x

Ch. 11 - Custom Functions

function [ outputArgs ] = MainFunctionName (inputArgs)

% statements go here

end

function [ outputArgs ] = LocalFunctionName (inputArgs)

% statements go here

end

- You are allowed to have nested functions in MATLAB

Anonymous Function:

- fctName = @(argList) expression
- Ex: RaisedCos = @(angle) (cosd(angle))^2;

- global variables - can be accessed from anywhere in the file
- Persistent variables
- persistent variable - only known to function where it was declared, maintains value between calls to function.

- Recursion - base case, decreasing case, ending case
- nargin - evaluates to the number of arguments the function was called with

Ch. 12 - Plotting

- plot(xArray, yArray)
- refer to help plot for more extensive documentation, chapter 12 only briefly covers plotting

plot - Line plot

yyaxis - Enables plotting with y-axes on both left and right sides

loglog - Line plot with logarithmic x and y axes

semilogy - Line plot with linear x and logarithmic y axes

semilogx - Line plot with logarithmic x and linear y axes

stairs - Stairstep graph

axis - sets the aspect ratio of x and y axes, ticks etc.

grid - adds a grid to the plot

gtext - allows positioning of text with the mouse

text - allows placing text at specified coordinates of the plot

xlabel labels the x-axis

ylabel labels the y-axis

title sets the graph title

figure(n) makes figure number n the current figure

hold on allows multiple plots to be superimposed on the same axes

hold off releases hold on current plot and allows a new graph to be drawn

close(n) closes figure number n

subplot(a, b, c) creates an a x b matrix of plots with c the current figure

orient specifies the orientation of a figure

Animated plot example:

for j = 1:360

pause(0.02)

plot3(x(j),y(j),z(j),'O')

axis([minx maxx miny maxy minz maxz]);

hold on;

end

Ch. 13 - Strings

stringArray = string(inputArray) - converts the input array inputArray to a string array

number = strLength(stringIn) - returns the number of characters in the input string

stringArray = strings(n, m) - returns an n-by-m array of strings with no characters,

- doing strings(sz) returns an array of strings with no characters, where sz defines the size.

charVec1 = char(stringScalar) char(stringScalar) converts the string scalar stringScalar into a character vector charVec1.

charVec2 = char(numericArray) char(numericArray) converts the numeric array numericArray into a character vector charVec2 corresponding to the Unicode transformation format 16 (UTF-16) code.

stringOut = string1 + string2 - combines the text in both strings

stringOut = join(textSrray) - consecutive elements of array joined, placing a space character between them.

stringOut = blanks(n) - creates a string of n blank characters

stringOut = strcar(string1, string2) - horizontally concatenates strings in arrays.

sprintf(formatSpec, number) - for printing out strings

- strExp = sprintf("%0.6e", pi)

stringArrayOur = compose(formatSpec, arrayIn) - formats data in arrayIn.

lower(string) - converts to lowercase

upper(string) - converts to uppercase

num2str(inputArray, precision) - returns a character vector with the maximum number of digits specified by precision

mat2str(inputMat, precision), converts matrix into a character vector.

numberOut = sscanf(inputText, format) - convert inputText according to format specifier

str2double(inputText)

str2num(inputChar)

strcmp(string1, string2)

strcmpi(string1, string2) - case-insensitive comparison

strncmp(str1, str2, n) - first n characters

strncmpi(str1, str2, n) - case-insensitive comparison of first n characters.

isstring(string) - logical output

isStringScalar(string) - logical output

ischar(inputVar) - logical output

iscellstr(inputVar) - logical output.

isstrprop(stringIn, categoryString) - returns a logical array of the same size as stringIn, with true at indices where elements of the stringIn belong to specified category:

iskeyword(stringIn) - returns true if string is a keyword in the matlab language

isletter(charVecIn)

isspace(charVecIn)

ischar(charVecIn)

contains(string1, testPattern) - boolean outputs if string contains a specific pattern

startsWith(string1, testPattern) - also logical output

endsWith(string1, testPattern) - also logical output

strfind(stringIn, pattern) - returns INDEX of each occurence in array

number = count(stringIn, patternSeek) - returns the number of occurences of string scalar in the string scalar stringIn.

strip(strArray) - removes all consecutive whitespace characters from begining and end of each string in Array, with side argument 'left', 'right', 'both'.

pad(stringIn) - pad with whitespace characters, can also specify where.

replace(stringIn, oldStr, newStr) - replaces all occurrences of oldStr in string array stringIn with newStr.

replaceBetween(strIn, startStr, endStr, newStr)

strrep(origStr, oldSubsr, newSubstr) - searches original string for substric, and if found, replace with new substric.

insertAfter(stringIn, startStr, newStr) - insert a new string afte the substring specified by startStr.

insertBefore(stringIn, endPos, newStr)

extractAfter(stringIn, startStr)

extractBefore(stringIn, startPos)

split(stringIn, delimiter) - divides string at whitespace characters.

splitlines(stringIn, delimiter)

It's frustrating when a long function or script runs and prints unexpected outputs to the command window. The line producing those outputs can be difficult to find.

Run this line of code before running the script or function. Execution will pause when the line is hit and the file will open to that line. Outputs that are intentionaly displayed by functions such as disp() or fprintf() will be ignored.

dbstop if unsuppressed output

To turn this off,

dbclear if unsuppressed output

Time to time I need to filll an existing array with NaNs using logical indexing. A trick I discover is using arithmetics rather than filling. It is quite faster in some circumtances

A=rand(10000);

b=A>0.5;

tic; A(b) = NaN; toc

tic; A = A + 0./~b; toc;

If you know trick for other value filling feel free to post.

function drawframe(f)

% Create a figure

figure;

hold on;

axis equal;

axis off;

% Draw the roads

rectangle('Position', [0, 0, 2, 30], 'FaceColor', [0.5 0.5 0.5]); % Left road

rectangle('Position', [2, 0, 2, 30], 'FaceColor', [0.5 0.5 0.5]); % Right road

% Draw the traffic light

trafficLightPole = rectangle('Position', [-1, 20, 1, 0.2], 'FaceColor', 'black'); % Pole

redLight = rectangle('Position', [0, 20, 0.5, 1], 'FaceColor', 'red'); % Red light

yellowLight = rectangle('Position', [0.5, 20, 0.5, 1], 'FaceColor', 'black'); % Yellow light

greenLight = rectangle('Position', [1, 20, 0.5, 1], 'FaceColor', 'black'); % Green light

carBody = rectangle('Position', [2.5, 2, 1, 4], 'Curvature', 0.2, 'FaceColor', 'red'); % Body

leftWheel = rectangle('Position', [2.5, 3.0, 0.2, 0.2], 'Curvature', [1, 1], 'FaceColor', 'black'); % Left wheel

rightWheel = rectangle('Position', [3.3, 3.0, 0.2, 0.2], 'Curvature', [1, 1], 'FaceColor', 'black'); % Right wheel

leftFrontWheel = rectangle('Position', [2.5, 5.0, 0.2, 0.2], 'Curvature', [1, 1], 'FaceColor', 'black'); % Left wheel

rightFrontWheel = rectangle('Position', [3.3, 5.0, 0.2, 0.2], 'Curvature', [1, 1], 'FaceColor', 'black'); % Right wheel

% Set limits

xlim([-1, 8]);

ylim([-1, 35]);

% Animation parameters

carSpeed = 0.5; % Speed of the car

carPosition = 2; % Initial car position

stopPosition = 15; % Position to stop at the traffic light

isStopped = false; % Car is not stopped initially

%Animation loop

for t = 1:100

% Update traffic light: Red for 40 frames, yellow for 10 frames Green for 40 frames

if t <= 40

% Red light on, yellow and green off

set(redLight, 'FaceColor', 'red');

set(yellowLight, 'FaceColor', 'black');

set(greenLight, 'FaceColor', 'black');

elseif t > 40 && t <= 50

% Change to green light

set(redLight, 'FaceColor', 'black');

set(yellowLight, 'FaceColor', 'yellow');

set(greenLight, 'FaceColor', 'black');

else

% Back to red light

set(redLight, 'FaceColor', 'black');

set(yellowLight, 'FaceColor', 'black');

set(greenLight, 'FaceColor', 'green');

isStopped = false; % Allow car to move

end

%Move the car

if ~isStopped

carPosition = carPosition + carSpeed; % Move forward

if carPosition < stopPosition

%do nothing

else

isStopped = true;

end

else

% Gradually stop the car when red

if carPosition > stopPosition

carPosition = carPosition + carSpeed*(1-t/50); % Move backward until it reaches the stop position

end

end

if carPosition >= 25

carPosition = 25;

end

% Update car position

% set(carBody, 'Position', [carPosition, 2, 1, 0.5]);

set(carBody, 'Position', [2.5, carPosition, 1, 4]);

%set(carWindow, 'Position', [carPosition + 0.2, 2.4, 0.6, 0.2]);

%set(leftWheel, 'Position', [carPosition, 1.5, 0.2, 0.2]);

set(leftWheel, 'Position', [2.5, carPosition+1, 0.2, 0.2]);

% set(rightWheel, 'Position', [carPosition + 0.8, 1.5, 0.2, 0.2]);

set(rightWheel, 'Position', [3.3, carPosition+1, 0.2, 0.2]);

set(leftFrontWheel, 'Position', [2.5, carPosition+3, 0.2, 0.2]);

set(rightFrontWheel, 'Position', [3.3, carPosition+3, 0.2, 0.2]);

% Pause to control animation speed

pause(0.01);

end

hold off;

In the spirit of warming up for this year's minihack contest, I'm uploading a walkthrough for how to design an airship using pure Matlab script. This is commented and uncondensed; half of the challenge for the minihacks is how minimize characters. But, maybe it will give people some ideas.

The actual airship design is from one of my favorite original NES games that I played when I was a kid - Little Nemo: The Dream Master. The design comes from the intro of the game when Nemo sees the Slumberland airship leave for Slumberland:

(Snip from a frame of the opening scene in Capcom's game Little Nemo: The Dream Master, showing the Slumberland airship).

I spent hours playing this game with my two sisters, when we were little. It's fun and tough, but the graphics sparked the imagination. On to the code walkthrough, beginning with the color palette: these four colors are the only colors used for the airship:

c1=cat(3,1,.7,.4); % Cream color

c2=cat(3,.7,.1,.3); % Magenta

c3=cat(3,0.7,.5,.1); % Gold

c4=cat(3,.5,.3,0); % bronze

We start with the airship carriage body. We want something rectangular but smoothed on the corners. To do this we are going to start with the separate derivatives of the x and y components, which can be expressed using separate blocks of only three levels: [1, 0, -1]. You could integrate to create a rectangle, but if we smooth the derivatives prior to integrating we will get rounded edges. This is done in the following code:

% Binary components for x & y vectors

z=zeros(1,30);

o=ones(1,100);

% X and y vectors

x=[z,o,z,-o];

y=[1+z,1-o,z-1,1-o];

% Smoother function (fourier / circular)

s=@(x)ifft(fft(x).*conj(fft(hann(45)'/22,260)));

% Integrator function with replication and smoothing to form mesh matrices

u=@(x)repmat(cumsum(s(x)),[30,1]);

% Construct x and y components of carriage with offsets

x3=u(x)-49.35;

y3=u(y)+6.35;

y3 = y3*1.25; % Make it a little fatter

% Add a z-component to make the full set of matrices for creating a 3D

% surface:

z3=linspace(0,1,30)'.*ones(1,260)*30;

z3(14,:)=z3(15,:); % These two lines are for adding platforms

z3(2,:)=z3(3,:); % to the carriage, later.

Plotting x, y, and the top row of the smoothed, integrated, and replicated matrices x3 and y3 gives the following:

We now have the x and y components for a 3D mesh of the carriage, let's make it more interesting by adding a color scheme including doors, and texture for the trim around the doors. Let's also add platforms beneath the doors for passengers to walk on.

Starting with the color values, let's make doors by convolving points in a color-matrix by a door shaped function.

m=0*z3; % Image matrix that will be overlayed on carriage surface

m(7,10:12:end)=1; % Door locations (lower deck)

m(21,10:12:end)=1; % Door locations (upper deck)

drs = ones(9, 5); % Door shape

m=1-conv2(m,ones(9,5),'same'); % Applying

To add the trim, we will convolve matrix "m" (the color matrix) with a square function, and look for values that lie between the extrema. We will use this to create a displacement function that bumps out the -x, and -y values of the carriage surface in intermediary polar coordinate format.

rm=conv2(m,ones(5)/25,'same'); % Smoothing the door function

rm(~m)=0; % Preserving only the region around the doors

rds=0*m; % Radial displacement function

rds(rm<1&rm>0)=1; % Preserving region around doors

rds(m==0)=0;

rds(13:14,:)=6; % Adding walkways

rds(1:2,:)=6;

% Apply radial displacement function

[th,rd]=cart2pol(x3,y3);

[x3T,y3T]=pol2cart(th,(rd+rds)*.89);

If we plot the color function (m) and radial displacement function (rds) we get the following:

In the upper plot you can see the doors, and in the bottom map you can see the walk way and door trim.

Next, we are going to add some flags draped from the bottom and top of the carriage. We are going to recycle the values in "z3" to do this, by multiplying that matrix with the absolute value of a sine-wave, squished a bit with the soft-clip erf() function.

We add a keel to the airship carriage using a canonical sphere turned on its side, again using the soft-clip erf() function to make it roughly rectangular in x and y, and multiplying with a vector that is half nan's to make the top half transparent.

At this point, since we are beginning the plotting of the ship, we also need to create our hgtransform objects. These allow us to move all of the components of the airship in unison, and also link objects with pivot points to the airship, such as the propeller.

% Now we need some flags extending around the top and bottom of the

% carriage. We can do this my multiplying the height function (z3) with the

% absolute value of a sine-wave, rounded with a compression function

% (erf() in this case);

g=-z3.*erf(abs(sin(linspace(0,40*pi,260))))/4; % Flags

% Also going to add a slight taper to the carriage... gives it a nice look

tp=linspace(1.05,1,30)';

% Finally, plotting. Plot the carriage with the color-map for the doors in

% the cream color, than the flags in magenta. Attach them both to transform

% objects for movement.

% Set up transform objects. 2 moving parts:

% 1) The airship itself and all sub-components

% 2) The propellor, which attaches to the airship and spins on its axis.

hold on;

P=hgtransform('Parent',gca); % Ship

S=hgtransform('Parent',P); % Prop

surf(x3T.*tp,y3T,z3,c1.*m,'Parent',P);

surf(x3,y3,g,c2.*rd./rd, 'Parent', P);

surf(x3,y3,g+31,c2.*rd./rd, 'Parent', P);

axis equal

% Now add the keel of the airship. Will use a canonical sphere and the

% erf() compression function to square off.

[x,y,z]=sphere(99);

mk=round(linspace(-1,1).^2+.3); % This function makes the top half of the sphere nan's for transparency.

surf(50*erf(1.4*z),15*erf(1.4*y),13*x.*mk./mk-1,.5*c2.*z./z, 'Parent', P);

% The carriage is done. Now we can make the blimp above it.

We haven't adjusted the shading of the image yet, but you can see the design features that have been created:

Next, we start working on the blimp. This is going to use a few more vertices & faces. We are going to use a tapered cylinder for this part, and will start by making the overlaid image, which will have 2 colors plus radial rings, circles, and squiggles for ornamentation.

M=525; % Blimp (matrix dimensions)

N=700;

% Assign the blimp the cream and magenta colors

t=122; % Transition point

b=ones(M,N,3); % Blimp color map template

bc=b.*c1; % Blimp color map

bc(:,t+1:end-t,:)=b(:,t+1:end-t,:).*c2;

% Add axial rings around blimp

l=[.17,.3,.31,.49];

l=round([l,1-fliplr(l)]*N); % Mirroring

lnw=ones(1,N); % Mask

lnw(l)=0;

lnw=rescale(conv(lnw,hann(7)','same'));

bc=bc.*lnw;

% Now add squiggles. We're going to do this by making an even function in

% the x-dimension (N, 725) added with a sinusoidal oscillation in the

% y-dimension (M, 500), then thresholding.

r=sin(linspace(0, 2*pi, M)*10)'+(linspace(-1, 1, N).^6-.18)*15;

q=abs(r)>.15;

r=sin(linspace(0, 2*pi, M)*12)'+(abs(linspace(-1, 1, N))-.25)*15;

q=q.*(abs(r)>.15);

% Now add the circles on the blimp. These will be spaced evenly in the

% polar angle dimension around the axis. We will have 9. To make the

% circles, we will create a cone function with a peak at the center of the

% circle, and use thresholding to create a ring of appropriate radius.

hs=[1,.75,.5,.25,0,-.25,-.5,-.75,-1]; % Axial spacing of rings

% Cone generation and ring loop

xy= @(h,s)meshgrid(linspace(-1, 1, N)+s*.53,(linspace(-1, 1, M)+h)*1.15);

w=@(x,y)sqrt(x.^2+y.^2);

for n=1:9

h=hs(n);

[xx,yy]=xy(h,-1);

r1=w(xx,yy);

[xx,yy]=xy(h,1);

r2=w(xx,yy);

b=@(x,y)abs(y-x)>.005;

q=q.*b(.1,r1).*b(.075,r1).*b(.1,r2).*b(.075,r2);

end

The figures below show the color scheme and mask used to apply the squiggles and circles generated in the code above:

Finally, for the colormap we are going to smooth the binary mask to avoid hard transitions, and use it to to add a "puffy" texture to the blimp shape. This will be done by diffusing the mask iteratively in a loop with a non-linear min() operator.

% 2D convolution function

ff=@(x)circshift(ifft2(fft2(x).*conj(fft2(hann(7)*hann(7)'/9,M,N))),[3,3]);

q=ff(q); % Smooth our mask function

hh=rgb2hsv(q.*bc); % Convert to hsv: we are going to use the value

% component for radial displacement offsets of the

% 3D blimp structure.

rd=hh(:,:,3); % Value component

for n=1:10

rd=min(rd,ff(rd)); % Diffusing the value component for a puffy look

end

rd=(rd+35)/36; % Make displacements effects small

% Now make 3D blimp manifold using "cylinder" canonical shape

[x,y,z]=cylinder(erf(sin(linspace(0,pi,N)).^.5)/4,M-1); % First argument is the blimp taper

[t,r]=cart2pol(x, y);

[x2,y2]=pol2cart(t, r.*rd'); % Applying radial displacment from mask

s=200;

% Plotting the blimp

surf(z'*s-s/2, y2'*s, x2'*s+s/3.9+15, q.*bc,'Parent',P);

Notice that the parent of the blimp surface plot is the same as the carriage (e.g. hgtransform object "P"). Plotting at this point using flat shading and adding some lighting gives the image below:

Next, we need to add a propeller so it can move. This will include the creation of a shaft using the cylinder() function. The rest of the pieces (the propeller blades, collars and shaft tip) all use the same canonical sphere with distortions applied using various math functions. Note that when the propeller is made it is linked to hgtransform object "S" rather than "P." This will allow the propeller to rotate, but still be joined to the airship.

% Next, the propeller. First, we start with the shaft. This is a simple

% cylinder. We add an offset variable and a scale variable to move our

% propeller components around, as well.

shx = -70; % This is our x-shifter for components

scl = 3; % Component size scaler

[x,y,z]=cylinder(1, 20); % Canonical cylinder for prop shaft.

p(1)=surf(-scl*(z-1)*7+shx,scl*x/2,scl*y/2,0*x+c4,'Parent',P); % Prop shaft

% Now the propeller. This is going to be made from a distorted sphere.

% The important thing here is that it is linked to the "S" hgtransform

% object, which will allow it to rotate.

[x,y,z]=sphere(50);

a=(-1:.04:1)';

x2=(x.*cos(a)-y/3.*sin(a)).*(abs(sin(a*2))*2+.1);

y2=(x.*sin(a)+y/3.*cos(a));

p(2)=surf(-scl*y2+shx,scl*x2,scl*z*6,0*x+c3,'Parent',S);

% Now for the prop-collars. You can see these on the shaft in the NES

% animation. These will just be made by using the canonical sphere and the

% erf() activation function to square it in the x-dimension.

g=erf(z*3)/3;

r=@(g)surf(-scl*g+shx,scl*x,scl*y,0*x+c3,'Parent',P);

r(g);

r(g-2.8);

r(g-3.7);

% Finally, the prop shaft tip. This will just be the sphere with a

% taper-function applied radially.

t=1.7*cos(0:.026:1.3)'.^2;

p(3)=surf(-(z*2+2)*scl + shx,x.*t*scl,y.*t*scl, 0*x+c4,'Parent',P);

Now for some final details including the ropes to the blimp, a flag hung on one of the ropes, and railings around the walkways so that passengers don't plummet to their doom. This will make use of the ad-hoc "ropeG" function, which takes a 3D vector of points and makes a conforming cylinder around it, so that you get lighting functions etc. that don't work on simple lines. This function is added to the script at the end to do this:

% Rope function for making a 3D curve have thickness, like a rope.

% Inputs:

% - xyz (3D curve vector, M points in 3 x M format)

% - N (Number of radial points in cylinder function around the curve

% - W (Width of the rope)

%

% Outputs:

% - xf, yf, zf (Matrices that can be used with surf())

function [xf, yf, zf] = RopeG(xyz, N, W)

% Canonical cylinder with N points in circumference

[xt,yt,zt] = cylinder(1, N);

% Extract just the first ring and make (W) wide

xyzt = [xt(1, :); yt(1, :); zt(1, :)]*W;

% Get local orientation vector between adjacent points in rope

dxyz = xyz(:, 2:end) - xyz(:, 1:end-1);

dxyz(:, end+1) = dxyz(:, end);

vcs = dxyz./vecnorm(dxyz);

% We need to orient circle so that its plane normal is parallel to

% xyzt. This is a kludgey way to do that.

vcs2 = [ones(2, size(vcs, 2)); -(vcs(1, :) + vcs(2, :))./(vcs(3, :)+0.01)];

vcs2 = vcs2./vecnorm(vcs2);

vcs3 = cross(vcs, vcs2);

p = @(x)permute(x, [1, 3, 2]);

rmats = [p(vcs3), p(vcs2), p(vcs)];

% Create surface

xyzF = pagemtimes(rmats, xyzt) + permute(xyz, [1, 3, 2]);

% Outputs for surf format

xf = squeeze(xyzF(1, :, :));

yf = squeeze(xyzF(2, :, :));

zf = squeeze(xyzF(3, :, :));

end

Using this function we can define the ropes and balconies. Note that the balconies simply recycle one of the rows of the original carriage surface, defining the outer rim of the walkway, but bumping up in the z-dimension.

cb=-sqrt(1-linspace(1, 0, 100).^2)';

c1v=[linspace(-67, -51)', 0*ones(100,1),cb*30+35];

c2v=[c1v(:,1),c1v(:,2),(linspace(1,0).^1.5-1)'*15+33];

c3v=c2v.*[-1,1,1];

[xr,yr,zr]=RopeG(c1v', 10, .5);

surf(xr,yr,zr,0*xr+c2,'Parent',P);

[xr,yr,zr]=RopeG(c2v', 10, .5);

surf(xr,yr,zr,0*zr+c2,'Parent',P);

[xr,yr,zr]=RopeG(c3v', 10, .5);

surf(xr,yr,zr,0*zr+c2,'Parent',P);

% Finally, balconies would add a nice touch to the carriage keep people

% from falling to their death at 10,000 feet.

[rx,ry,rz]=RopeG([x3T(14, :); y3T(14,:); 0*x3T(14,:)+18]*1.01, 10, 1);

surf(rx,ry,rz,0*rz+cat(3,0.7,.5,.1),'Parent',P);

surf(rx,ry,rz-13,0*rz+cat(3,0.7,.5,.1),'Parent',P);

And, very last, we are going to add a flag attached to the outer cable. Let's make it flap in the wind. To make it we will recycle the z3 matrix again, but taper it based on its x-value. Then we will sinusoidally oscillate the flag in the y dimension as a function of x, constraining the y-position to be zero where it meets the cable. Lastly, we will displace it quadratically in the x-dimension as a function of z so that it lines up nicely with the rope. The phase of the sine-function is modified in the animation loop to give it a flapping motion.

h=linspace(0,1);

sc=10;

[fx,fz]=meshgrid(h,h-.5);

F=surf(sc*2.5*fx-90-2*(fz+.5).^2,sc*.3*erf(3*(1-h)).*sin(10*fx+n/5),sc*fz.*h+25,0*fx+c3,'Parent',P);

Plotting just the cables and flag shows:

Putting all the pieces together reveals the full airship:

A note about lighting: lighting and material properties really change the feel of the image you create. The above picture is rendered in a cartoony style by setting the specular exponent to a very low value (1), and adding lots of diffuse and ambient reflectivity as well. A light below the airship was also added, albeit with lower strength. Settings used for this plot were:

shading flat

view([0, 0]);

L=light;

L.Color = [1,1,1]/4;

light('position', [0, 0.5, 1], 'color', [1,1,1]);

light('position', [0, 1, -1], 'color', [1, 1, 1]/5);

material([1, 1, .7, 1])

set(gcf, 'color', 'k');

axis equal off

What about all the rest of the stuff (clouds, moon, atmospheric haze etc.) These were all (mostly) recycled bits from previous minihack entries. The clouds were made using power-law noise as explained in Adam Danz' blog post. The moon was borrowed from moonrun, but with an increased number of points. Atmospheric haze was recycled from Matlon5. The rest is just camera angles, hgtransform matrix updates, and updating alpha-maps or vertex coordinates.

Finally, the use of hann() adds the signal processing toolbox as a dependency. To avoid this use the following anonymous function:

hann = @(x)-cospi(linspace(0,2,x)')/2+.5;

Create a struct arrays where each struct has field names "a," "b," and "c," which store different types of data. What efficient methods do you have to assign values from individual variables "a," "b," and "c" to each struct element? Here are five methods I've provided, listed in order of decreasing efficiency. What do you think?

Create an array of 10,000 structures, each containing each of the elements corresponding to the a,b,c variables.

num = 10000;

a = (1:num)';

b = string(a);

c = rand(3,3,num);

Here are the methods;

%% method1

t1 =tic;

s = struct("a",[], ...

"b",[], ...

"c",[]);

s1 = repmat(s,num,1);

for i = 1:num

s1(i).a = a(i);

s1(i).b = b(i);

s1(i).c = c(:,:,i);

end

t1 = toc(t1);

%% method2

t2 =tic;

for i = num:-1:1

s2(i).a = a(i);

s2(i).b = b(i);

s2(i).c = c(:,:,i);

end

t2 = toc(t2);

%% method3

t3 =tic;

for i = 1:num

s3(i).a = a(i);

s3(i).b = b(i);

s3(i).c = c(:,:,i);

end

t3 = toc(t3);

%% method4

t4 =tic;

ct = permute(c,[3,2,1]);

t = table(a,b,ct);

s4 = table2struct(t);

t4 = toc(t4);

%% method5

t5 =tic;

s5 = struct("a",num2cell(a),...

"b",num2cell(b),...

"c",squeeze(mat2cell(c,3,3,ones(num,1))));

t5 = toc(t5);

%% plot

bar([t1,t2,t3,t4,t5])

xtickformat('method %g')

ylabel("time(second)")

yline(mean([t1,t2,t3,t4,t5]))

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.

D.R. Kaprekar was a self taught recreational mathematician, perhaps known mostly for some numbers that bear his name.

Today, I'll focus on Kaprekar's constant (as opposed to Kaprekar numbers.)

The idea is a simple one, embodied in these 5 steps.

1. Take any 4 digit integer, reduce to its decimal digits.

2. Sort the digits in decreasing order.

3. Flip the sequence of those digits, then recompose the two sets of sorted digits into 4 digit numbers. If there were any 0 digits, they will become leading zeros on the smaller number. In this case, a leading zero is acceptable to consider a number as a 4 digit integer.

4. Subtract the two numbers, smaller from the larger. The result will always have no more than 4 decimal digits. If it is less than 1000, then presume there are leading zero digits.

5. If necessary, repeat the above operation, until the result converges to a stable result, or until you see a cycle.

Since this process is deterministic, and must always result in a new 4 digit integer, it must either terminate at either an absorbing state, or in a cycle.

For example, consider the number 6174.

7641 - 1467

We get 6174 directly back. That seems rather surprising to me. But even more interesting is you will find all 4 digit numbers (excluding the pure rep-digit nmbers) will always terminate at 6174, after at most a few steps. For example, if we start with 1234

4321 - 1234

8730 - 0378

8532 - 2358

and we see that after 3 iterations of this process, we end at 6174. Similarly, if we start with 9998, it too maps to 6174 after 5 iterations.

9998 ==> 999 ==> 8991 ==> 8082 ==> 8532 ==> 6174

Why should that happen? That is, why should 6174 always drop out in the end? Clearly, since this is a deterministic proces which always produces another 4 digit integer (Assuming we treat integers with a leading zero as 4 digit integers), we must either end in some cycle, or we must end at some absorbing state. But for all (non-pure rep-digit) starting points to end at the same place, it seems just a bit surprising.

I always like to start a problem by working on a simpler problem, and see if it gives me some intuition about the process. I'll do the same thing here, but with a pair of two digit numbers. There are 100 possible two digit numbers, since we must treat all one digit numbers as having a "tens" digit of 0.

N = (0:99)';

Next, form the Kaprekar mapping for 2 digit numbers. This is easier than you may think, since we can do it in a very few lines of code on all possible inputs.

Ndig = dec2base(N,10,2) - '0';

Nmap = sort(Ndig,2,'descend')*[10;1] - sort(Ndig,2,'ascend')*[10;1];

I'll turn it into a graph, so we can visualize what happens. It also gives me an excuse to employ a very pretty set of tools in MATLAB.

G2 = graph(N+1,Nmap+1,[],cellstr(dec2base(N,10,2)));

plot(G2)

Do you see what happens? All of the rep-digit numbers, like 11, 44, 55, etc., all map directly to 0, and they stay there, since 0 also maps into 0. We can see that in the star on the lower right.

G2cycles = cyclebasis(G2)

G2cycles{1}

All other numbers eventually end up in the cycle:

G2cycles{2}

That is

81 ==> 63 ==> 27 ==> 45 ==> 09 ==> and back to 81

looping forever.

Another way of trying to visualize what happens with 2 digit numbers is to use symbolics. Thus, if we assume any 2 digit number can be written as 10*T+U, where I'll assume T>=U, since we always sort the digits first

syms T U

(10*T + U) - (10*U+T)

So after one iteration for 2 digit numbers, the result maps ALWAYS to a new 2 digit number that is divisible by 9. And there are only 10 such 2 digit numbers that are divisible by 9. So the 2-digit case must resolve itself rather quickly.

What happens when we move to 3 digit numbers? Note that for any 3 digit number abc (without loss of generality, assume a >= b >= c) it almost looks like it reduces to the 2 digit probem, aince we have abc - cba. The middle digit will always cancel itself in the subtraction operation. Does that mean we should expect a cycle at the end, as happens with 2 digit numbers? A simple modification to our previous code will tell us the answer.

N = (0:999)';

Ndig = dec2base(N,10,3) - '0';

Nmap = sort(Ndig,2,'descend')*[100;10;1] - sort(Ndig,2,'ascend')*[100;10;1];

G3 = graph(N+1,Nmap+1,[],cellstr(dec2base(N,10,2)));

plot(G3)

This one is more difficult to visualize, since there are 1000 nodes in the graph. However, we can clearly see two disjoint groups.

We can use cyclebasis to tell us the complete story again.

G3cycles = cyclebasis(G3)

G3cycles{:}

And we see that all 3 digit numbers must either terminate at 000, or 495. For example, if we start with 181, we would see:

811 - 118

963 - 369

954 - 459

It will terminate there, forever trapped at 495. And cyclebasis tells us there are no other cycles besides the boring one at 000.

What is the maximum length of any such path to get to 495?

D3 = distances(G3,496) % Remember, MATLAB uses an index origin of 1

D3(isinf(D3)) = -inf; % some nodes can never reach 495, so they have an infinite distance

plot(D3)

The maximum number of steps to get to 495 is 6 steps.

find(D3 == 6) - 1

So the 3 digit number 100 required 6 iterations to eventually reach 495.

shortestpath(G3,101,496) - 1

I think I've rather exhausted the 3 digit case. It is time now to move to the 4 digit problem, but we've already done all the hard work. The same scheme will apply to compute a graph. And the graph theory tools do all the hard work for us.

N = (0:9999)';

Ndig = dec2base(N,10,4) - '0';

Nmap = sort(Ndig,2,'descend')*[1000;100;10;1] - sort(Ndig,2,'ascend')*[1000;100;10;1];

G4 = graph(N+1,Nmap+1,[],cellstr(dec2base(N,10,2)));

plot(G4)

cyclebasis(G4)

ans{:}

And here we see the behavior, with one stable final point, 6174 as the only non-zero ending state. There are no circular cycles as we had for the 2-digit case.

How many iterations were necessary at most before termination?

D4 = distances(G4,6175);

D4(isinf(D4)) = -inf;

plot(D4)

The plot tells the story here. The maximum number of iterations before termination is 7 for the 4 digit case.

find(D4 == 7,1,'last') - 1

shortestpath(G4,9986,6175) - 1

Can you go further? Are there 5 or 6 digit Kaprekar constants? Sadly, I have read that for more than 4 digits, things break down a bit, there is no 5 digit (or higher) Kaprekar constant.

We can verify that fact, at least for 5 digit numbers.

N = (0:99999)';

Ndig = dec2base(N,10,5) - '0';

Nmap = sort(Ndig,2,'descend')*[10000;1000;100;10;1] - sort(Ndig,2,'ascend')*[10000;1000;100;10;1];

G5 = graph(N+1,Nmap+1,[],cellstr(dec2base(N,10,2)));

plot(G5)

cyclebasis(G5)

ans{:}

The result here are 4 disjoint cycles. Of course the rep-digit cycle must always be on its own, but the other three cycles are also fully disjoint, and are of respective length 2, 4, and 4.

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;

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

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

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