Trouble with the "double" numeric data type
5 次查看(过去 30 天)
显示 更早的评论
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 个评论
Dyuman Joshi
2023-5-31
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))
IntN = IntN1-IntN2;
TN(s,r) = IntN;
disp(TN)
回答(1 个)
VBBV
2023-5-31
编辑:VBBV
2023-5-31
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;
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!