How to simplify infinite double matrix summation in Matlab.

I need to calculate the following sum:

I used the following code:

beta = 1.5; A = [1 0;0 2]; R = [0 3; 1 1];
syms k l t
f = ((-A).^k./factorial(k)).*((factorial(k+l).*(-R).^l)/gamma((2-beta).*l...
    + 2.*k + 2)).*((t.^((2-beta).*l+2.*k+1))/factorial(l));
F = vpasum(vpasum(f,k,0,inf),l,0,inf);
G = subs(F,t,1)

G = 

This yields G as a symbolic sum but I need it as numeric sum. I tried it using symsum function but that too yields a symbolic answer. Using, subs function for t=1 won't change it into a numeric value. Any help would be highly appriciated.

[Torsten]

All parameters of your double sum are known. So why do you work with symbolic variables if you want numeric results ?

[Walter Roberson]

@Torsten

... because the summations are infinite.

[Torsten]

... because the summations are infinite.

What does it help if no analytical expression for the symbolic sums exists ?

[Walter Roberson]

(-R).^l

I suspect that is

(-R)^l

and

(-A).^k

I suspect that is

(-A)^k

[Hilal Ahmad Bhat]

@Walter Roberson Yes, you are right but that won't make any difference in this case.

Answer 1

[Walter Roberson]

beta = 1.5; A = [1 0;0 2]; R = [0 3; 1 1];
syms k l t
f = ((-A)^k./factorial(k)).*((factorial(k+l).*(-R)^l)/gamma((2-beta).*l...
    + 2.*k + 2)).*((t.^((2-beta).*l+2.*k+1))/factorial(l));
F = symsum(symsum(f,k,0,inf),l,0,inf);
G = subs(F,t,1)

G = 

vpa(G)

ans = 

[Torsten]

I think there should also be matrix multiplication instead of elementwise multiplication between (-A)^k and (-R)^l, shouldn't it ?

[Walter Roberson]

Good point.

beta = 1.5; A = [1 0;0 2]; R = [0 3; 1 1];
syms k l t
f = ((-A)^k./factorial(k))*((factorial(k+l)*(-R)^l)/gamma((2-beta).*l...
    + 2.*k + 2)).*((t.^((2-beta).*l+2.*k+1))/factorial(l));
F = symsum(symsum(f,k,0,inf),l,0,inf);
G = subs(F,t,1)

G = 

vpa(G)

ans = 

[Hilal Ahmad Bhat]

@Walter Roberson @Torsten

Thank you very much to both of you.

[Hilal Ahmad Bhat]

I was trying to solve the question using while loops. I wrote the following code:

clear all; beta = 1.5; A = [1 0;0 2]; R = [0 3; 1 1];  %#ok<CLALL>
inn_total = zeros(2); inn_total(:) = inf;
out_total = zeros(2);
l=0; k=0; t=1;

while sqrt(trace(inn_total*inn_total')) > eps % Norm of inn_total should get as much small as it pleases for the convergence of outer summation.
    inn_total = zeros(2);
    term=zeros(2); term(:) = inf;
    while sqrt(trace(term*term')) > eps % Norm of term should get as much small as it pleases for the convergence of inner summation.
        term = ((-A)^k/factorial(k))*((factorial(k+l)*(-R)^l)/...
               gamma((2-beta)*l+ 2*k + 2))*((t^((2-beta)*l+2*k+1))/...
               factorial(l));
        inn_total = inn_total + term;
        l = l+1;
    end    
    out_total = out_total + inn_total;
    k = k+1;
end
out_total %#ok<NOPTS>

out_total = 2×2

3.1688 -3.6799 -1.2266 1.9422

[k,l] %#ok<NOPTS>

ans = 1×2

2 71

The idea that I am using in the convergence of the series is mathematical, i.e., converges if norm of each term after some finite 'n' must be less than some predefined small value (eps in my case). For inner sum, while loop stops after value of 'Lth term' is less than 'eps' and the outer while loop stops when the norm of inner sum gets less than 'eps'.

But the values, I get from this method, differ by approx 0.6 than the following two methods - one using the symsum and vpa (@Walter Roberson's method), and the other using the loops (@Torsten's method).

% Walter Roberson Method
beta = 1.5; A = [1 0;0 2]; R = [0 3; 1 1];
syms k l t
f1 = ((-A)^k/factorial(k))*((factorial(k+l)*(-R)^l)/gamma((2-beta)*l...
    + 2*k + 2))*((t^((2-beta)*l+2*k+1))/factorial(l));
F = symsum(symsum(f1,k,0,inf),l,0,inf);
G = subs(F,t,1);
vpa(G,5)

ans = 

% Torsten Method
beta = 1.5; A = [1 0;0 2]; R = [0 3; 1 1]; t = 1;
N = 4;
M = 36;
mat = zeros(2);
for s1 = 0:N
    for s2 = 0:M
        mat = mat + f(s1,s2,beta,A,R,t);
    end
end
mat

mat = 2×2

2.5141 -2.8682 -0.7202 1.2210

function mat = f(k,l,beta,A,R,t)
  mat = (-A)^k/factorial(k)*(-R)^l/factorial(l)*factorial(k+l)...
        *t^((2-beta)*l+2*k+1)/gamma((2-beta)*l+2*k+2);
end

Why does the values differ when I use while loop? Am I missing anything or doing something extra in my code. I need this to run in while loop because it is very useful sometimes. Thanks in advance.

[Torsten]

beta = 1.5; A = [1 0;0 2]; R = [0 3; 1 1]; t = 1;
mat = zeros(2);
mat_inner = Inf(2);
tol = 1e-6;
s1 = 0;
while norm(mat_inner,'fro') > tol
    mat_inner = zeros(2);
    for s2 = 0:s1
        mat_inner = mat_inner + f(s1-s2,s2,beta,A,R,t);
    end
    mat = mat + mat_inner;
    s1 = s1 + 1;
end
s1

s1 = 45

format long
mat

mat = 2x2

2.514129941450880 -2.868263680288298 -0.720173713227138 1.221004006453150

function mat = f(k,l,beta,A,R,t)
  mat = (-1)^(k+l)*A^k/factorial(k)*R^l/factorial(l)*factorial(k+l)...
        *t^((2-beta)*l+2*k+1)/gamma((2-beta)*l+2*k+2);
end

[Hilal Ahmad Bhat]

@Torsten

This is awesome. Thank you very much.

Answer 2

[Torsten]

Are you sure that the matrix potentials for A and R are meant elementwise ?

beta = 1.5; A = [1 0;0 2]; R = [0 3; 1 1]; t = 1;
N = 50;
M = N;
mat = zeros(2);
for s1 = 0:N
    for s2 = 0:M
        mat = mat + f(s1,s2,beta,A,R,t);
    end
end
mat

mat = 2×2

0.8415 0.2856 0.5560 0.4253

function mat = f(k,l,beta,A,R,t)
  mat = (-A).^k/factorial(k).*(-R).^l/factorial(l)*factorial(k+l)...
        *t^((2-beta)*l+2*k+1)/gamma((2-beta)*l+2*k+2);
end

[Hilal Ahmad Bhat]

I need infinite sum. Using loops, we are able to calculate the sum upto finite terms only.

[Hilal Ahmad Bhat]

Oh! It is my mistake. The matrix potential is not elementwise. It is the meant to be the matrix multiplication (i.e., rows × colums). Whatever be the case, I am stuck at finding the infinite sum. I think while loop might help but I don't know how to do it with double summation.

[Torsten]

I doubt you will find an analytical expression for your double sum - that's what summing up to infinity would mean. Thus you will have to compute the sum numerically.

Infinite sums are numerically evaluated by building finite sums up to an index N until the result changes only neclectably if N is increased.

That's what I did.

If you meant usual matrix multiplication, replace

function mat = f(k,l,beta,A,R,t)
  mat = (-A).^k/factorial(k).*(-R).^l/factorial(l)*factorial(k+l)...
        *t^((2-beta)*l+2*k+1)/gamma((2-beta)*l+2*k+2);
end

by

function mat = f(k,l,beta,A,R,t)
  mat = (-A)^k/factorial(k)*(-R)^l/factorial(l)*factorial(k+l)...
        *t^((2-beta)*l+2*k+1)/gamma((2-beta)*l+2*k+2);
end