Trouble with the "double" numeric data type

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?

[Dyuman Joshi]

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]

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

Answer 1

[VBBV]

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;

[FRANCESCO FABRIS]

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.