# Trouble with the "double" numeric data type

3 views (last 30 days)
FRANCESCO FABRIS on 31 May 2023
Commented: FRANCESCO FABRIS on 31 May 2023
I have the following code - that is extracted from a "for r = 1:length(mu)" cycle not working properly - in which we compute the results for r = 24:
clear;
syms x t
r=24;
deltamu = 0.5;
sigma = 5;
lambda = 1;
Cstar = 7;
mu = [-4:deltamu:10];
P = zeros(1,length(mu));
TP = zeros(1,length(mu));
TN = zeros(1,length(mu));
s=1;
sigmae = 5;
exp_t = exp(-t.^2);
tmax = +Inf;
tmin = ((mu(r)+lambda.*sigma.^2-x)./(sqrt(2).*sigma));
erfc = (2./sqrt(pi)).*int(exp_t,t,tmin,tmax);
f_lambda = (lambda./2).*exp((lambda./2).*(2.*mu(r)+lambda.*sigma.^2-2.*x)).*erfc;
P(s,r) = double(int(f_lambda,x,Cstar,+Inf));
phi_e = 1./(sigmae.*sqrt(2.*pi)).*exp(-(1/2).*((x-Cstar)./sigmae).^2);
PHI = int(phi_e,x,-Inf,x);
fxP = f_lambda.* PHI;
TP(s,r) = double(int(fxP,x,Cstar,+Inf));
fxN1 = f_lambda;
fxN2 = f_lambda.*PHI;
IntN1 = double(int(fxN1,x,-Inf,Cstar));
IntN2 = double(int(fxN2,x,-Inf,Cstar));
IntN = IntN1-IntN2;
TN(s,r) = double(IntN)
It works well, computing the 24th component of the vector TN
TN =
Columns 1 through 20
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Columns 21 through 29
0 0 0 0.2809 0 0 0 0 0
Nevertheless, quite surprisingly, if we set instead r = 25, we obtain the following message:
"Error using symengine
Unable to convert expression containing remaining symbolic function calls into double array. Argument must be expression that evaluates to number. Error in sym/double (line 729)
Xstr = mupadmex('symobj::double', S.s, 0);
Related documentation"
This error pops up for all values 25 <= r <=29, where 29 is its maximum allowed.
At a closer inspection, it appears that the problem is associated with the last "double" calculation, i.e.
IntN2 = double(int(fxN2,x,-Inf,Cstar));
that is done correctly for r=24, but gives an error with r >=25. These are the two values obtained:
****************
r=24
>> int(fxN2,x,-Inf,Cstar)
ans =
int((146084674377193476124496039018555*2^(1/2)*pi*exp(20 - x)*(erf((2^(1/2)*(x - 65/2))/10) + 1)*(erf(2^(1/2)*(x/10 - 7/10)) + 1))/2596148429267413814265248164610048, x, -Inf, 7)
double(int(fxN2,x,-Inf,Cstar))
ans =
0.1042
****************
r=25
>> int(fxN2,x,-Inf,Cstar)
ans =
int((146084674377193476124496039018555*2^(1/2)*pi*exp(41/2 - x)*(erf((2^(1/2)*(x - 33))/10) + 1)*(erf(2^(1/2)*(x/10 - 7/10)) + 1))/2596148429267413814265248164610048, x, -Inf, 7)
double(int(fxN2,x,-Inf,Cstar))
Error using symengine
Unable to convert expression containing remaining symbolic function calls into double array. Argument must be expression that evaluates to number.
Error in sym/double (line 729)
Xstr = mupadmex('symobj::double', S.s, 0);
Related documentation
****************
It is evident that the only difference between the two version of the function fxN2 - that will be integrated - is exp(20 - x) vs exp(41/2 - x) and (x - 65/2) vs (x - 33) inside the erf function.
Where is the trick?
##### 2 CommentsShow 1 older commentHide 1 older comment
FRANCESCO FABRIS on 31 May 2023
Thank you so much, Dyuma, for your nice suggestion. It works fine and solves my problem! It remains the mistery of why the transaction r=24 -> 25 prevent the correct calculation

VBBV on 31 May 2023
Edited: VBBV on 31 May 2023
The below version of code runs without errors.
clear;
syms x t
% r=24;
deltamu = 0.5;
sigma = 5;
lambda = 1;
Cstar = 7;
mu = [-4:deltamu:10]
P = zeros(1,length(mu));
TP = zeros(1,length(mu));
TN = zeros(1,length(mu));
s=1;
sigmae = 5;
exp_t = exp(-t.^2);
tmax = +Inf;
for r = 1:length(mu)
tmin = ((mu(r)+lambda.*(sigma.^2)-x)./(sqrt(2).*sigma));
erfc = (2./sqrt(pi)).*int(exp_t,t,tmin,tmax);
f_lambda = (lambda./2).*exp((lambda./2).*(2.*mu(r)+lambda.*(sigma.^2)-x)).*erfc;
P(s,r) = double(int(f_lambda,x,Cstar,+Inf));
phi_e = 1./(sigmae.*sqrt(2.*pi)).*exp(-(1/2).*((x-Cstar)./sigmae).^2);
PHI = int(phi_e,x,-Inf,x);
fxP = f_lambda.* PHI;
TP(s,r) = double(int(fxP,x,Cstar,+Inf));
fxN1 = f_lambda;
fxN2 = f_lambda.*PHI;
IntN1 = double(int(fxN1,x,-Inf,Cstar));
IntN2 = double(int(fxN2,x,-Inf,Cstar));
IntN = IntN1-IntN2;
TN(s,r) = double(IntN)
end
if I understand correctly, do you mean writing this line in code as
%---------------------->>------->>
tmin = ((mu(r)+lambda.*(sigma.^2)-x)./(sqrt(2).*sigma));
%------------------------------------------------------------>>------->>
f_lambda = (lambda./2).*exp((lambda./2).*(2.*mu(r)+lambda.*(sigma.^2)-x)).*erfc;
##### 1 CommentShow NoneHide None
FRANCESCO FABRIS on 31 May 2023
Dear VBBV, thank you for your attempt. Your code works fine, but you are using a different f_lambda.
Mine is
f_lambda = (lambda./2).*exp((lambda./2).*(2.*mu(r)+lambda.*sigma.^2-2.*x)).*erfc;
while yours is
f_lambda = (lambda./2).*exp((lambda./2).*(2.*mu(r)+lambda.*(sigma.^2)-x)).*erfc;
It is another mystery why this little difference generates a logical error. Thank you in any case.

### Categories

Find more on Operators and Elementary Operations in Help Center and File Exchange

R2023a

### Community Treasure Hunt

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

Start Hunting!