"Unable to compute number of steps from 0 to n by 1" error for symbolic integration

Im trying to do numerical integration with the following function that have several parameters. My issue is that whenever I assign n as symbolic symbol, I got the error "Unable to compute number of steps from 0 to n by 1". I Also tried using symsum instead of arrayfunc and I got the same error.

n=60;
syms u_f positive 
syms lam_f positive 
syms z positive 
syms alpha positive 
syms sigma positive
syms ep positive
syms z positive
syms b positive
syms d positive
assume (ep <1 & ep>0)
%assume(z < u_f )
%syms n 
% assume (n, {'postive', 'integer'})

nCk = @(n,kVec,z)arrayfun(@(k)lam_f*exp(-lam_f*(alpha-b)).*(nchoosek(n,k)*(1*ep).^(k).*((1-ep)*(0.5+0.5*erf((z-u_f-d)/(sqrt(2)*sigma)))).^(n-(k)).*(0.5+0.5*erf((z-alpha-u_f-d)/(sqrt(2)*sigma)))^k),kVec);
%F(z,u_f,lam_f,sigma,ep)=int(nCk(n,0:n,z),alpha,0, z);
F(z,u_f,lam_f,sigma,ep,b,d)=vpaintegral(nCk(n,0:n,z),alpha,b, z);
result=vpaintegral( 1-sum(vpa(F(z,0.4,7.1111,0.009,0.1,2,2))),0,Inf,'ArrayValued',true,'RelTol', 1e-5, 'AbsTol', 1e-4)

My end goal is writing a user defined function like the following:

function result=norm_expo_gang_function_mean_handle(n,uu,lamm,oo,epp,bb,dd)
%(u,lam,o,b,d,ep,sz1,sz2)
syms F(z,u_f,lam_f,sigma_o,b,d,ep)
syms u_f positive 
syms lam_f positive 
syms z positive 
syms alpha_o positive 
syms sigma_o positive
syms ep positive
syms z positive
syms b positive
syms d positive
assume (ep <1 & ep>0)
assume(z,'positive')
assume(lam_f,'positive');
assume(u_f,'positive');
assume(sigma_o,'positive');
assume(b,'positive');
assume(d,'positive');
%assume(n,{'positive','integer'});

assume (ep <1 & ep>0)
nCk = @(n,kVec,z)arrayfun(@(k)lam_f*exp(-lam_f*(alpha-b)).*(nchoosek(n,k)*(1*ep).^(k).*((1-ep)*(0.5+0.5*erf((z-u_f-d)/(sqrt(2)*sigma)))).^(n-(k)).*(0.5+0.5*erf((z-alpha-u_f-d)/(sqrt(2)*sigma)))^k),kVec);
%nCk = @(n,kVec,z)arrayfun(@(k)lam_f*exp(-lam_f*(alpha-b)).*(nchoosek(n,k)*(1*ep).^(k).*((1-ep)*(0.5+0.5*erf((z-u_f-d)/(sqrt(2)*sigma)))).^(n-(k)).*(0.5+0.5*erf((z-alpha-u_f-d)/(sqrt(2)*sigma)))^k),kVec);
%F(z,u_f,lam_f,sigma,ep)=int(nCk(n,0:n,z),alpha,0, z);

F(z,u_f,lam_f,sigma_o,ep,b,d)=vpaintegral(nCk(n,0:n,z),alpha,b, z);
result=vpaintegral( 1-sum(vpa(F(z,uu,lamm,oo,epp,bb,dd,n))),0,Inf,'ArrayValued',true,'RelTol', 1e-5, 'AbsTol', 1e-4);
%F(z,u_f,lam_f,sigma_o,b,d,ep,n) = @(n,kVec,z)arrayfun(@(k)lam_f*exp(-lam_f*(alpha_o-b)).*(nchoosek(n,k)*(1*ep).^(k).*((1-ep)*(0.5+0.5*erf((z-u_f-d)/(sqrt(2)*sigma_o)))).^(n-(k)).*(0.5+0.5*erf((z-alpha_o-u_f-d)/(sqrt(2)*sigma_o)))^k),kVec);

end

Answer 1

[Walter Roberson]

syms z positive
syms b positive
%...
F(z,u_f,lam_f,sigma,ep,b,d)=vpaintegral(nCk(n,0:n,z),alpha,b, z);

You are trying to vpaintegral() with a symbolic lower and upper bound. When you vpaintegral() the bounds must be definite numbers

[Walter Roberson]

Q = @(v) sym(v);
n=60;
syms u_f positive 
syms lam_f positive 
syms z positive 
syms alpha positive 
syms sigma positive
syms ep positive
syms z positive
syms b positive
syms d positive
assume (ep <1 & ep>0)
%assume(z < u_f )
%syms n 
% assume (n, {'postive', 'integer'})

nCk = @(n,kVec,z)arrayfun(@(k)lam_f*exp(-lam_f*(alpha-b)).*(nchoosek(n,k)*(1*ep).^(k).*((1-ep)*(0.5+0.5*erf((z-u_f-d)/(sqrt(2)*sigma)))).^(n-(k)).*(0.5+0.5*erf((z-alpha-u_f-d)/(sqrt(2)*sigma)))^k),kVec);
%F(z,u_f,lam_f,sigma,ep)=int(nCk(n,0:n,z),alpha,0, z);
F(z,u_f,lam_f,sigma,ep,b,d)=int(nCk(n,Q(0:n),z),alpha,b, z, 'hold', true)
result=int( 1-sum((F(z,Q(0.4),Q(7.1111),Q(0.009),Q(0.1),Q(2),Q(2)))),0,Inf,'hold', true)
%release(result)

the release() step times out here on MATLAB Answers.

[Muna Tageldin]

I get an errror "Too many arguments". I'm not sure what caused this.

[Walter Roberson]

What step are you encountering the error on?

My guess at the moment is that you are using a MATLAB version that is too old to support the 'hold' option ?

[Muna Tageldin]

Yes, the error was on the hold operation and I use MATLAB 2020.

[Walter Roberson]

Q = @(v) sym(v);
n=60;
syms u_f positive 
syms lam_f positive 
syms z positive 
syms alpha positive 
syms sigma positive
syms ep positive
syms z positive
syms b positive
syms d positive
assume (ep <1 & ep>0)
%assume(z < u_f )
%syms n 
% assume (n, {'postive', 'integer'})

nCk = @(n,kVec,z)arrayfun(@(k)lam_f*exp(-lam_f*(alpha-b)).*(nchoosek(n,k)*(1*ep).^(k).*((1-ep)*(0.5+0.5*erf((z-u_f-d)/(sqrt(2)*sigma)))).^(n-(k)).*(0.5+0.5*erf((z-alpha-u_f-d)/(sqrt(2)*sigma)))^k),kVec);
F(z,u_f,lam_f,sigma,ep,b,d)=int(nCk(n,Q(0:n),z),alpha,b, z)
result=int( 1-sum((F(z,Q(0.4),Q(7.1111),Q(0.009),Q(0.1),Q(2),Q(2)))),0,Inf)

This will take a fair while -- roughly 5 minutes.

After that you may wish to

R = vpa(result)

[Muna Tageldin]

Why did you use int in your code? Vpaintegral is suited for numerical integrals (no analytical solution possible) based on this page https://www.mathworks.com/help/symbolic/sym.vpaintegral.html

I'm confused on the use of the whole code.

[Walter Roberson]

nested vpaintegral calls are less accurate than vpa of a complete integral.

Also you were attempting to vpaintegral with symbolic boundaries, which is not permitted

[Muna Tageldin]

Ok. I understand your point. I have a question, how would you display the numerical value of the result. I tried to do vpa(ans) but all I get is the unevaluated integral

int((9*exp(7111/500 - (7111*z)/1000))/200 - (9*erf(2^(1/2)*((500*z)/9 - 2200/9)))/200 - (9*erf((400*2^(1/2))/3)*exp(7111/500 - (7111*z)/1000))/200 + (9*exp(62580895872001/2000000000000 - (7111*z)/1000)*erf((500000^(1/2)*z)/9 - (4400575991*500000^(1/2))/9000000000))/200 - (63999*erf((500*2^(1/2)*(z - 12/5))/9)*((1000*erf(2^(1/2)*((500*z)/9 - 2200/9)))/7111 + (1000*erf((400*2^(1/2))/3)*exp(7111/500 - (7111*z)/1000))/7111 - (1000*exp(62580895872001/2000000000000 - (7111*z)/1000)*erf((500000^(1/2)*z)/9 - (4400575991*500000^(1/2))/9000000000))/7111 - (1000*erf((800191997*500000^(1/2))/3000000000)*exp(62580895872001/2000000000000 - (7111*z)/1000))/7111))/200000 + (9*erf((800191997*500000^(1/2))/3000000000)*exp(62580895872001/2000000000000 - (7111*z)/1000))/200 + (63999*erf((500*2^(1/2)*(z - 12/5))/9)*((1000*exp(7111/500 - (7111*z)/1000))/7111 - 1000/7111))/200000 - int((7111*exp(7111/500 - (7111*alpha)/1000)*(erf((500*2^(1/2)*(alpha - z + 12/5))/9)/2 - 1/2)^2)/100000, alpha, 2, z) + (575991*((1000*exp(-(7111*z)/1000)*exp(7111/500))/7111 - 1000/7111)*(erf((500*2^(1/2)*(z - 12/5))/9)/2 + 1/2)^2)/100000 + 191/200, z, 0, Inf)