Are there alternative ways to write codes to calculate a_p ?

Question

Rohitashya 2024-10-23，14:45

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2161695-are-there-alternative-ways-to-write-codes-to-calculate-a_p

评论： Torsten 2024-10-23，17:56

在 MATLAB Online 中打开

Here

given

m= 0:n+1

n is the degree of Bspline which can be any value

p = 0:n
k=k0:k1;
k0 =ceil(d-((n+1)/2));
k1 =floor(d+((n+1)/2));
num_points = 1; % Change this if you need more points
% Generate random variable d uniformly in the range (-0.5, 0.5)
d = -0.5 + (0.5 - (-0.5)) * rand(num_points, 1);

mu is the unit step function or heaviside function

I tried to develop the function a_P

function a_p = calculate_ap(n, d)
% n: Degree of the polynomial
% d: Interpolation point as a random variable with |d| < 0.5
% Define the range for p (0 to n)
p_values = 0:n;
% Calculate the lower and upper bounds of k
k0 = ceil(d - (n + 1) / 2);
k1 = floor(d + (n + 1) / 2);
% Define the range for k
k_values = k0:k1;
% Initialize a matrix to store a_p(k) values
a_p = zeros(length(p_values), length(k_values));
% Loop through each value of p to compute a_p for each k
for p_idx = 1:length(p_values)
    p = p_values(p_idx);
    
    % Calculate the binomial coefficient (n choose p)
    binomial_np = nchoosek(n, p);
    
    % Calculate the factorial term 1/n!
    n_factorial_term = 1 / factorial(n);
    
    % Loop through each value of k
    for k_idx = 1:length(k_values)
        k = k_values(k_idx);
        
        % Initialize the sum for this combination of p and k
        sum_result = 0;
        
        % Sum over m from 0 to n+1
        for m = 0:n+1
            % Calculate the binomial coefficient (n+1 choose m)
            binomial_n1m = nchoosek(n+1, m);
            
            % Calculate the term (-1)^m
            sign_term = (-1)^m;
            
            % Calculate the term (n+1)/2 - m
            difference_term = (n+1) / 2 - m;
            
            % Calculate the term (difference_term)^(n-p)
            power_term = difference_term^(n - p);
            
            % Calculate the Heaviside function mu(p - m + (n+1)/2)
            mu_term = heaviside(p - m + (n+1)/2);
            
            % Add the current term to the sum
            sum_result = sum_result + binomial_n1m * sign_term * power_term * mu_term;
        end
        
        % Calculate a_p(k) for the current p and k
        a_p(p_idx, k_idx) = n_factorial_term * binomial_np * sum_result;
    end
end

For solving this function

n= input("Enter the degree:"); %Degree
p=0:n;
a_p = calculate_ap(n,p);
disp(a_p)

solution is :

for n= 3

0.6667 0.6667 0.6667 0.6667 0.6667

-1.0000 -1.0000 -1.0000 -1.0000 -1.0000

0.5000 0.5000 0.5000 0.5000 0.5000

0 0 0 0 0

>> Is this function written correctly does it need any change ? Please could you suggest?

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Torsten 2024-10-23，17:56

I wonder why a_p depends on k, but the variable k does not appear in its definition.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Subhajyoti 2024-10-23，15:58

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2161695-are-there-alternative-ways-to-write-codes-to-calculate-a_p#answer_1536385

在 MATLAB Online 中打开

Hi @Rohitashya,

You can try precomputing values and vectorized operations to efficiently calculate the given equation. Additionally, dynamic programming paradigm helps in avoiding repeated calculations and optimizes code efficiency.

Here, in the following implementation, I have utilized precomputed factorial and binomials to compute the equation.

MAX_N = 10000;
% preallocate memory for factorial values upto MAX_N
factorial_values = cumprod(1:MAX_N);
function y = fact(n)
    global MAX_N, factorial_values;
    if n==0
        y = 1;
    else 
        y = factorial_values(n);
    end
end

% preallocate 2-dimentions array to precompute binomial coefficients
binomial_values = zeros(MAX_N, MAX_N);
% precompute binomial coefficients using dynamic programming
for n = 1:MAX_N
    for k = 0:n
        if k == 0 || k == n
            binomial_values(n+1, k+1) = 1;
        else
            binomial_values(n+1, k+1) = binomial_values(n, k) + binomial_values(n, k+1);
        end
    end
end
function nCr(n, r)
    global binomial_values;
    binomial_values(n+1, r+1)
end

function a_p = calculate_ap(n, p)
    ans = 0;
    for m = 0:n+2
         sign_term = (-1)^m;
    
         base = ((n+1)/2 - m);
         power_term = (base)^(n-p);
    
         % Calculate the Heaviside function mu(p - m + (n+1)/2)
         mu = heaviside(p - m + (n+1)/2);
    
         ans = ans + nCr(n+1, m) * sign_term * power_term * mu;
    end
   ans = ans / (fact(n-p) * fact(p));
end

For repetitive calculations, you can try to further optimize the code by memorizing the “heaviside” function and utilising binary exponentiation algorithm.

Additionally, you can refer to the following resources to know more about using common techniques used to improve code run-time:

3 个评论
显示 1更早的评论隐藏 1更早的评论

Rohitashya 2024-10-23，16:11

MAX_N = 10000;

% preallocate memory for factorial values upto MAX_N

factorial_values = cumprod(1:MAX_N);

function y = fact(n)

global MAX_N, factorial_values;

if n==0

y = 1;

else

y = factorial_values(n);

end

end could you explain this part a bit ?

Rohitashya 2024-10-23，16:14

Please explain in details about the frst two parts

请先登录，再进行评论。

Are there alternative ways to write codes to calculate a_p ?

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

回答（1 个）

3 个评论
显示 1更早的评论隐藏 1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Are there alternative ways to write codes to calculate a_p ?

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

回答（1 个）

3 个评论 显示 1更早的评论隐藏 1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

3 个评论
显示 1更早的评论隐藏 1更早的评论