Clear Filters
Clear Filters

Trouble with the "double" numeric data type

2 views (last 30 days)
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 Comments
Dyuman Joshi
Dyuman Joshi on 31 May 2023
Not all integrals can be evaluated by MATLAB symbolically (or other symbolic solvers for that matter). In such cases you can use vpaintegral to numerically approximate the value of the integral -
syms x t
r=25;
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(vpaintegral(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(vpaintegral(fxP,x,Cstar,+Inf));
fxN1 = f_lambda;
fxN2 = f_lambda.*PHI;
IntN1 = double(vpaintegral(fxN1,x,-Inf,Cstar));
IntN2 = double(vpaintegral(fxN2,x,-Inf,Cstar))
IntN2 = 0.1042
IntN = IntN1-IntN2;
TN(s,r) = IntN;
disp(TN)
Columns 1 through 19 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Columns 20 through 29 0 0 0 0 0.2809 0 0 0 0 0
FRANCESCO FABRIS
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

Sign in to comment.

Answers (1)

VBBV
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 Comment
FRANCESCO FABRIS
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.

Sign in to comment.

Categories

Find more on Creating and Concatenating Matrices in Help Center and File Exchange

Tags

Products


Release

R2023a

Community Treasure Hunt

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

Start Hunting!