How to solve following equation in matlab

2 views (last 30 days)
Syed
Syed on 12 Feb 2018
Commented: Walter Roberson on 13 Feb 2018
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)
  2 Comments
Syed
Syed on 13 Feb 2018
@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...

Sign in to comment.

Answers (1)

Walter Roberson
Walter Roberson on 13 Feb 2018
f6_35 = double(sum(subs(f5,jj,0:35)))
gives as precise an answer in double precision as going to 200 does.
  5 Comments
Walter Roberson
Walter Roberson on 13 Feb 2018
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)))

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!