Evaluating numerical value of symbolic integral

I need numerical value of Eb_avg, instead matlab is giving symbolic expression and is not performing the integrals

Following is my code:

clc; 
clear all; 
close all;
neta=0.75;                               
rho=0.5;  
T=1e-3;
P_No_dB=15;                                                      
P_No_end=size(P_No_dB);
in=P_No_end(1,2);
P=10.^(P_No_dB./10);       
sigma_sq=1;                        
Ebmax=100*1e-3; 
d0=1;                              
m=3;                              
var_hsd=(d0^(-m));  
d1=2;                           
m=3;                              
var_hse=(d1^(-m));  
d2=2;                          
m=3;                             
var_hed=(d2^(-m));                                                              
R=2;                             
SNRth=2^R;                    
count=0;
beta=(2^R-1).*sigma_sq./(P.*(1-rho));
lambda_0=1/(var_hsd);
lambda_1=1/(var_hse);
lambda_2=1/(var_hed);
Eb=[];
for j=1:in                                                     
syms x y    
fun=(T*lambda_2.*((x-y.*(1-rho))./(y.*(1-rho))).*sigma_sq.*( expint(lambda_2.*(T.*(x-y.*(1-rho)).*sigma_sq)./(y.*(1-rho).*(T.*neta.*rho.*P.*y +Ebmax))) -expint(lambda_2.*((x-y.*(1-rho)).*sigma_sq)./(y.*(1-rho).*(neta.*rho.*P.*y))) ) -T.*neta.*rho.*P.*y.*( exp(-lambda_2.*(T.*(x-y.*(1-rho)).*sigma_sq)./(y.*(1-rho).*(T.*neta.*rho.*P.*y +Ebmax))) - exp(-lambda_2.*((x-y.*(1-rho)).*sigma_sq)./(y.*(1-rho).*(neta.*rho.*P.*y))) )).*lambda_1.*lambda_0.*exp(-lambda_1.*y).*exp(-lambda_0.*x);
fun1=int(fun,y,beta,x/(1-rho));
Eb_avg=int(fun1,x,0,Inf);
Eb=[Eb Eb_avg];
count=count+1
end

[Torsten]

Use the numerical integrator "integral2".

[madhan ravi]

In addition to Torsten's comment aren't you calculating the same thing over and over again? I haven't paid enough attention to details.

[Walter Roberson]

for j=1:in would normally be doing the same calculations over and over again, but the upper bound happens to be defined as the second dimension of a scalar, so in = 1

Answer 1

[Walter Roberson]

change both int() into vpaintegral().

You have a function in two variables, and it does not look to have a closed form in either of the two variables individually, so the first int() is going to return unresolved. It looks like MATLAB is not able to find a closed form for the 2D integral. So switch to numeric integration on formulae by using vpaintegral()

Note: make sure you replace both int() calls.

... It still won't be fast. Integrating to infinity seldom is.

[Amit Patel]

Problem still persists using vpaintegral(), matlab is returning symbolic expression.

I have variable upper limit in the first integral (fun1) , so how can I use integral2 for numerical integration

[Amit Patel]

if I use integral to accommodate the variable upper limit-

fun=@(x,y) (T*lambda_2.*((x-y.*(1-rho))./(y.*(1-rho))).*sigma_sq.*( expint(lambda_2.*(T.*(x-y.*(1-rho)).*sigma_sq)./(y.*(1-rho).*(T.*neta.*rho.*P.*y +Ebmax))) -expint(lambda_2.*((x-y.*(1-rho)).*sigma_sq)./(y.*(1-rho).*(neta.*rho.*P.*y))) ) -T.*neta.*rho.*P.*y.*( exp(-lambda_2.*(T.*(x-y.*(1-rho)).*sigma_sq)./(y.*(1-rho).*(T.*neta.*rho.*P.*y +Ebmax))) - exp(-lambda_2.*((x-y.*(1-rho)).*sigma_sq)./(y.*(1-rho).*(neta.*rho.*P.*y))) )).*lambda_1.*lambda_0.*exp(-lambda_1.*y).*exp(-lambda_0.*x);
  fun1=@(x) integral(fun,beta,x/(1-rho),'ArrayValued',1);
  Eb_avg= integral(fun1,0,Inf,'ArrayValued',1);
  

Then i am getting following error-

  Not enough input arguments.

Error in
AvgEb_Consumption_25Jan>@(x,y)(T*lambda_2.*((x-y.*(1-rho))./(y.*(1-rho))).*sigma_sq.*(expint(lambda_2.*(T.*(x-y.*(1-rho)).*sigma_sq)./(y.*(1-rho).*(T.*neta.*rho.*P.*y+Ebmax)))-expint(lambda_2.*((x-y.*(1-rho)).*sigma_sq)./(y.*(1-rho).*(neta.*rho.*P.*y))))-T.*neta.*rho.*P.*y.*(exp(-lambda_2.*(T.*(x-y.*(1-rho)).*sigma_sq)./(y.*(1-rho).*(T.*neta.*rho.*P.*y+Ebmax)))-exp(-lambda_2.*((x-y.*(1-rho)).*sigma_sq)./(y.*(1-rho).*(neta.*rho.*P.*y))))).*lambda_1.*lambda_0.*exp(-lambda_1.*y).*exp(-lambda_0.*x)

Error in integralCalc/iterateArrayValued (line 156)
                fxj = FUN(t(1)).*w(1);

Error in integralCalc/vadapt (line 130)
            [q,errbnd] = iterateArrayValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
        [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 88)
Q = integralCalc(fun,a,b,opstruct);

Error in AvgEb_Consumption_25Jan>@(x)integral(fun,beta,x/(1-rho),'ArrayValued',1)

Error in integralCalc/iterateArrayValued (line 156)
                fxj = FUN(t(1)).*w(1);

Error in integralCalc/vadapt (line 130)
            [q,errbnd] = iterateArrayValued(u,tinterval,pathlen);

Error in integralCalc (line 83)
        [q,errbnd] = vadapt(@AToInfInvTransform,interval);

Error in integral (line 88)
Q = integralCalc(fun,a,b,opstruct);

Error in AvgEb_Consumption_25Jan (line 80)
  Eb_avg= integral(fun1,0,Inf,'ArrayValued',1);
  

What am I missing?

[Torsten]

fun = @(x,y)...;  % all parameters in "fun" have to be scalars !!!!
ymax = @(x)x/(1-rho);
Eb_avg = integral2(fun,0,Inf,beta,ymax);

[Walter Roberson]

you cannot use integral with aa function if two variables .

when you use integral2 then the second set of bounds can be function handle of function of one variable .

integral2(fun, 0, inf, beta, @(xx) xx/(1-rho) )

[Amit Patel]

Thanks , Integral issue is resolved but integral is not converging.

[Walter Roberson]

One thing I notice is that the integral of y is 3/50*sqrt(10) lower bound, 2*x upper bound, but with x starting from 0, 2*x can be less than 3/50*sqrt(10) so you are doing "backwards" integration for that stretch. I recommend checking that bound as that hints that possibly you are crossing a discontinuity.

[Amit Patel]

Yes , thanks for pointing out that issue of bound, I have resolved that and now integral is converging.