How to solve following equation in matlab

1 次查看（过去 30 天）
Syed
Syed 2018-2-12
评论： Walter Roberson 2018-2-13
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 个评论
Walter Roberson
Walter Roberson 2018-2-12
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
Syed 2018-2-13
@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...

请先登录，再进行评论。

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

回答（1 个）

Walter Roberson
Walter Roberson 2018-2-13
f6_35 = double(sum(subs(f5,jj,0:35)))
gives as precise an answer in double precision as going to 200 does.
  5 个评论
显示 3更早的评论
Syed
Syed 2018-2-13
@Walter Roberson, I am using R2015b About Gamma function, I am able to use it on other occasions in Matlab all fine.
Walter Roberson
Walter Roberson 2018-2-13
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)))

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Common Operations 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by