How to solve following equation in matlab

Hi all, please see the attached figure which is my equation I am trying to solve using Matlab.

I have the code using syms function (see below). I believe f6 should yield final result for me but instead it shows something of the following kind. I used functions like 'vpa' and 'simplify' but still cannot yield answer. See f6 first and then the code

f6 = (3884161996073403*symsum((5^jj*symsum((-1)^(jj - m)/factorial(jj - m), m, 0, jj))/factorial(jj), jj, 0, Inf))/576460752303423488

The Code is as follows

   clc;
  clear all;
  close all;

Pc = 45-30;
Pc_lin = 10^(Pc/10);

Po = 3-30;
Po_lin = 10^(Po/10);

lambda_m = 6e-6
alpha_l = 2;
alpha_o = 4;
varepsilon_2 = 0.5;
Xd = 5;

% T_mj = e^(-K)*K^j*(1/(factorial(j) * factorial(j-m)) )

beta = (2*pi/alpha_o)/sin(2*pi/alpha_o);

syms jj m

K = 5; 
    T = -5;
    T_lin = 10^(T/10);
    Cd = pi*(T_lin*varepsilon_2*Xd^alpha_l)^(2/alpha_o)*(lambda_m*(Pc_lin/Po_lin)^(2/alpha_o))*beta;
          A = Cd*T_lin^(2/alpha_o);
    B = exp(-1*Cd*T_lin^(2/alpha_o));         
    % T1
    syms jj m z ii
    f1 = (-1)^(jj-m)/(factorial(z-ii)*factorial(ii)) * gamma(2*ii/alpha_o)/gamma(2*ii/alpha_o-(jj-m)+1)
    f2 = symsum(f1, ii, 1, z)
    f3 = symsum(f2,z,1,jj-m);
    f4 = (-1)^(jj-m)*exp(-K)*K^jj/(factorial(jj)*factorial(jj-m))
    f5 = symsum(f4, m, 0, jj)
    f6 = symsum(f5, jj,0,inf)

[Walter Roberson]

You are not necessarily doing anything wrong. A lot of symbolic summations do not have known closed form solutions. I cannot read your image clearly, but with what I see, I speculate that you might not be able to do better than a solution that replaces one of the summations with a hypergeom call, and I am not certain that is possible.

[Syed]

@Walter Roberson, Any thing to solve this equation using symb? I tried running the forloop for this equation from 1 to 150 and it took 40 minutes and still the first run was not yet complete. I had to abondon that (since there were nearly 10 runs in total).

I tried the same thing with symsum using limit from 1 to 200 (since when j > 150, 1/j! is near to Zero. But then I got the following error (please see the image).

I tried to lower down the limit from 150 to 10 to check if this 'Singularity' can be resolved, but yet again, the same error...

Answer 1

[Walter Roberson]

f6_35 = double(sum(subs(f5,jj,0:35)))

gives as precise an answer in double precision as going to 200 does.

[Syed]

:( received the same singularity message again (after a long long wait)

[Walter Roberson]

Which MATLAB release are you using? I tested in R2017b. Note that I left out the assignment to f6, instead doing the f6_35 finite summation. Response is fairly quick.

[Walter Roberson]

The infinite summation for f6 does not appear to have a closed form representation. It comes out as a constant times

sum(5^jj*GAMMA(jj+1, -1)/GAMMA(jj+1)^2, jj = 0 .. infinity)

where GAMMA with two arguments is the incomplete GAMMA function.

[Syed]

@Walter Roberson, I am using R2015b About Gamma function, I am able to use it on other occasions in Matlab all fine.

[Walter Roberson]

It works for me when I test in R2015b. For clarity I am testing with

Pc = 45-30;
Pc_lin = 10^(Pc/10);

Po = 3-30;
Po_lin = 10^(Po/10);

lambda_m = 6e-6
alpha_l = 2;
alpha_o = 4;
varepsilon_2 = 0.5;
Xd = 5;

% T_mj = e^(-K)*K^j*(1/(factorial(j) * factorial(j-m)) )

beta = (2*pi/alpha_o)/sin(2*pi/alpha_o);

syms jj m

K = 5; 
    T = -5;
    T_lin = 10^(T/10);
    Cd = pi*(T_lin*varepsilon_2*Xd^alpha_l)^(2/alpha_o)*(lambda_m*(Pc_lin/Po_lin)^(2/alpha_o))*beta;
          A = Cd*T_lin^(2/alpha_o);
    B = exp(-1*Cd*T_lin^(2/alpha_o));         
    % T1
    syms jj m z ii
    f1 = (-1)^(jj-m)/(factorial(z-ii)*factorial(ii)) * gamma(2*ii/alpha_o)/gamma(2*ii/alpha_o-(jj-m)+1)
    f2 = symsum(f1, ii, 1, z)
    f3 = symsum(f2,z,1,jj-m);
    f4 = (-1)^(jj-m)*exp(-K)*K^jj/(factorial(jj)*factorial(jj-m))
    f5 = symsum(f4, m, 0, jj)

      f6_35 = double(sum(subs(f5,jj,0:35)))