How to convert symbolic output to numeric value in Symbolic Math Toolbox?

322 views (last 30 days)
Dear MATLAB users,
I tried to calculate following problem using Symbolic Math Toolbox.
Unfortunately, I have faced some problems using Symbolic Math Toolbox.
My code is as follows.
My problem is when i run this code, when the mu is small number, the P0, Pmu an Pc is calculated as numerical value.
But, starts from some points, the results are shown like 'round((513*vpasum(Inf/gamma(n + 181), n, 0, Inf))/2500)'.
When i tried to convert this kind of output to numerical value, it fails.
Therefore, is there any way to convert this kind of output to numerical value?
Thanks in advance!
syms n mu P0 Pmu Pc obj sigma lambda k
lambda=20
sigma=0.1
for mu=1:40
an=symprod(lambda/(mu+(k*sigma)),k,1,n);
P0=1/(1+symsum(an,n,1,Inf));
Pmu=mu*P0_record*vpa(symsum(((lambda^n)*gamma(1+(mu/sigma)))/((sigma^n)*(mu+n*sigma)*gamma(1+n+(mu/sigma))),n,0,Inf));
Pc=(1-P0)*(mu/(mu+sigma));
obj=Pmu*Pc;
P0_record(mu)=round(P0,4);
Pmu_record(mu)=round(Pmu,4);
Pc_record(mu)=round(Pc,4);
obj_record(mu)=round(obj,4);
end

Accepted Answer

Paul
Paul on 7 Feb 2023
Edited: Paul on 7 Feb 2023
My first inclination was to stay symbolic all the way through and then compute numerical values at the end.
syms lambda mu sigma positive
syms n k integer
a_n = symprod(lambda/(mu + k*sigma),k,1,n);
a_n = simplify(a_n,10)
a_n = 
P_0 = 1/(1 + symsum(a_n,n,1,inf));
P_0 = simplify(P_0,10)
P_0 = 
P_n = P_0*a_n
P_n = 
P_mu = symsum(P_n*(mu/(mu + n*sigma)),n,0,inf)
P_mu = 
P_c = (1 - P_0)*(mu/(mu + sigma))
P_c = 
Now evaluate using subs. For example
muval = (1:40).';
Pc_record = subs(P_c, [lambda sigma mu] , {20, 0.1 , muval});
Pmu_record = subs(P_mu, [lambda sigma mu] , {20, 0.1 , muval});
P0_record = subs(P_0, [lambda sigma mu] , {20, 0.1 , muval});
VPA to see the values. Need to increase digits to deal with very large terms in the expressions.
digits(100);
figure
h1 = subplot(211);hold on;plot(muval,real(vpa(Pc_record))),ylabel('real(Pc)')
h2 = subplot(212);hold on;plot(muval,imag(vpa(Pc_record))),ylabel('imag(Pc)');xlabel('mu')
figure
h3 = subplot(211);hold on;plot(muval,real(vpa(Pmu_record))),ylabel('real(Pmu)')
h4 = subplot(212);hold on;plot(muval,imag(vpa(Pmu_record))),ylabel('imag(Pmu)');xlabel('mu')
figure
h5 = subplot(211);hold on;plot(muval,real(vpa(P0_record))),ylabel('real(P0)')
h6 = subplot(212);hold on;plot(muval,imag(vpa(P0_record))),ylabel('imag(P0)'),xlabel('mu')
It looks like for large vaues of mu some error is creeping into the vpa calculations yielding a very, very, small imaginary component even with digits(100). I assume that imaginary part can be safely ignored. Or maybe use more digits? Or maybe the expressions can be further simplified symbolically before subbing in specific values. Sometimes that helps.
Alternatively, setting lambda and sigma ahead of the loop seems to work better, but slower. Couldn't run for all values of mu because it takes too long (takes much longer to run here than on my local machine)
digits(32) % reset to default
clear P0_record Pmu_record Pc_record
syms lambda mu sigma positive
syms n k integer
lambda = sym(20);
sigma = sym(0.1);
Can only use a few values of muval to stay under the Answers runtime limit
muval = [1:8:40 40];
for ii = 1:numel(muval)
mu = sym(muval(ii));
a_n = symprod(lambda/(mu + k*sigma),k,1,n);
a_n = simplify(a_n,10);
P_0 = 1/(1 + symsum(a_n,n,1,inf));
P_0 = simplify(P_0);
P_n = P_0*a_n;
P_mu = symsum(P_n*(mu/(mu + n*sigma)),n,0,inf);
P_c = (1 - P_0)*(mu/(mu + sigma));
Pc_record(ii) = double(P_c);
Pmu_record(ii) = double(P_mu);
P0_record(ii) = double(P_0);
end
plot(h1,muval,real(Pc_record),'-x'),grid
plot(h2,muval,imag(Pc_record),'-x'),grid
plot(h3,muval,real(Pmu_record),'-x'),grid
plot(h4,muval,imag(Pmu_record),'-x'),grid
plot(h5,muval,real(P0_record),'-x'),grid
plot(h6,muval,imag(P0_record),'-x'),grid

More Answers (2)

KSSV
KSSV on 7 Feb 2023
Read about double, vpasolve.

Amal Raj
Amal Raj on 7 Feb 2023
Hi,
You can use the vpa function to convert the symbolic expression to a numerical value with a specified number of digits.
Then go on with the calculation.
P0=1/(1+symsum(an,n,1,Inf));
Pmu=mu*P0_record*vpa(symsum(((lambda^n)*gamma(1+(mu/sigma)))/((sigma^n)*(mu+n*sigma)*gamma(1+n+(mu/sigma))),n,0,Inf));
Pc=(1-P0)*(mu/(mu+sigma));
% 32 significat digits
P0=vpa(1/(1+symsum(an,n,1,Inf)),32);
Pmu=vpa(mu*P0_record*vpa(symsum(((lambda^n)*gamma(1+(mu/sigma)))/((sigma^n)*(mu+n*sigma)*gamma(1+n+(mu/sigma))),n,0,Inf)),32);
Pc=vpa((1-P0)*(mu/(mu+sigma)),32);

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!