Representation of this vector T_d1 and Td_2 in MATLAB

For this information I tried to write the Td_2 in MATLAB

% Function to compute Stirling numbers of the first kind
function S = stirling1(n, k)
    if n == k
        S = 1;
    elseif k == 0
        S = 0;
    elseif n == 0
        S = 0;
    else
        S = stirling1(n-1, k-1) - (n-1) * stirling1(n-1, k);
    end
end
% Pre-allocate the matrix T_d^2
T_d2 = zeros(M, M);
% Fill the matrix according to the Stirling number formula
for i = 1:M
    for j = 1:i
        T_d2(i, j) = stirling1(i, j);  % Stirling number of the first kind
    end
end
disp('Matrix T_d^2:');
disp(T_d2);

The results show Matrix T_d^2:

1 0 0 0

-1 1 0 0

2 -3 1 0

-6 11 -6 1

Td_1 as % Pre-allocate the matrix T_d^1

T_d1 = zeros(M, M);
% Fill the matrix according to the formula
for i = 1:M
    for j = 1:i  % j ranges from 1 to i
        T_d1(i, j) = nchoosek(i-1, j-1) * ((-M + 1) / 2)^(i-j);
    end
end
disp('Matrix T_d^1:');
disp(T_d1);

The results are Matrix T_d^1:

1.0000 0 0 0

-1.5000 1.0000 0 0

2.2500 -3.0000 1.0000 0

-3.3750 6.7500 -4.5000 1.0000

Is it a correct method of calculation ? Please could you give suggestions or comments ?

Answer 1

[Ashok]

Hey @Rohitashya,

The provided code indeed calculates the and matrices based on the given formula. Here's an alternative approach that doesn't require declaring 'stirling1' as a separate function.

function Td1 = compute_Td1(M)
    % Create a matrix of binomial coefficients
    [I, J] = ndgrid(1:M, 1:M);
    binomials = zeros(M, M);
    for i = 1:M
        binomials(i, 1:i) = arrayfun(@(ii, jj) nchoosek(ii-1, jj-1), I(i, 1:i), J(i, 1:i)); 
    end
    % Create a matrix for powers
    powers = ((- (M - 1) / 2) .^ (I - J)) .* (J <= I);
    % Element-wise multiplication to get Td1
    Td1 = binomials .* powers;
end

function Td2 = compute_Td2(M)
    % Initialize Td2 with identity
    Td2 = eye(M+1);
    for n = 2:M+1
        for k = 2:n-1
            Td2(n, k) = Td2(n-1, k-1) - (n-2) * Td2(n-1, k); 
        end 
    end
    Td2 = Td2(2:M+1, 2:M+1);
end

% Example usage:
Td1 = compute_Td1(M);
disp('T_d^1:');
disp(Td1);

Td2 = compute_Td2(M);
disp('T_d^2:');
disp(Td2);

Additionally, for implementing cross products, consider using MATLAB's 'cross' function instead of nested for loops. You can find the relevant documentation here:

https://www.mathworks.com/help/releases/R2024b/matlab/ref/cross.html

I believe this will be of use!

[Rohitashya]

Thank You Ashok for suggesting this solution