Search code examples
matlabintegral

numerical integration of MATLAB on my function give me many numbers not a single number


Numerical integration of MATLAB on my function give me many numbers not a single number. I'm not good at MatLaB, so I don't know the reason. Is there anyone who can help me ?

this is my function.

function [ y ] = SucP( l,a,P,D,r,th );

d=2./a;
Ehd=integral(@(x)x.^d.*exp(-x),0,inf);
gam=gamma(1-d);
C1=1./(1+th.*(r/D).^a);
C2=1./(1+th.*P.*(r/D).^a);
hypgeo1=hypergeom([1,2],[2-d],th.*(r./D).^a.*C1);
hypgeo2=hypergeom([1,2],[2-d],th.*P.*(r./D).^a.*C2);
y=exp(-l.*pi.*(th).^d.*r.^2.*P.^d.*Ehd.*gam-C1+C2+th./(1-d).*r.^2.*(r./D).^(a-2).*(C1.^2.*hypgeo1-P.^d.*C2.^2.*hypgeo2));

end

And I want to integrate

fun=@(l,a,P,D,r,z)SucP(l,a,P,D,r,2.^z-1);
y=integral(@(z)fun(l,a,P,D,r,z),0,inf);

But this integration give me following result.

Columns 1 through 6

   0.999869167524854   0.998589370430984   0.994817933624792   0.987704371328770   0.976976845412355   0.962748430805626

  Columns 7 through 12

   0.945394762911001   0.925627034788835   0.904532141053543   0.883378353636190   0.863363119639727   0.845589211136565 ....

Columns 145 through 150

             NaN                 NaN                 NaN                 NaN                 NaN                 NaN


ans =

   NaN

Is there anyone who can help me ?


Solution

  • Try something along the lines of :

    1) You must get rid of those NaN values. The reason why they appear is the following :

    You're trying to integrate your function

    SucP(l,a,P,D,r,2.^z-1)
    

    From z=0 to infinity. But, from a certain value of z , 2^z-1 will be bigger than the biggest number matlab can deal with.

    What you need to do to avoid that is to change the bounds of your integration so that 2^z-1 never goes above realmax (which is the biggest value MATLAB can handle).

    The upper bound of your integration will be :

    Upperbound=log(realmax)/log(2)-4;
    

    2) You must change your declaration of the fun dunction you declare to :

    fun=@(z)SucP(l,a,P,D,r,2.^z-1);
    

    And then call :

    Y=integral(fun,0,log(realmax)/log(2)-4);