MATLAB dispalays Error in integralCa​lc/finalIn​putCheck

This is my code. I would like to plot p_c and p_d on the y-axis given the formula in the picture.

    big_X0=10^(-27.96/10); 
    r=5;
    P_tx_h = big_X0/(r^(-2));
    alpha = 4;
    delta = 2/alpha;
    P_c = 1;% 
    P_d = 1;
    T = 0:2:20;
    big_X=10.^((0:2:20)/10);
    r_dist = big_X/P_tx_h;
    lambda_c = 0.01;
    lambda_d = 0.01;
    mu=0.7854;
    Tplot= 10.^(0/10);
    %mu = (delta/(1-delta)).*Tplot.* hypergeom([1, 1-delta], 2-delta,-Tplot);%use if you have symbolic maths toolbox
    v = (Tplot.^delta).*gamma(1-delta).*gamma(1+delta);
    lambda_c_eq_intf = lambda_c + (lambda_d.*(P_d/P_c).^delta).*(v./mu);
    lambda_d_eq_intf = lambda_d + (lambda_c.*(P_c/P_d).^delta);
    f_r_dist=  @(r_dist) 2.*pi.*lambda_c.*r_dist.*exp(-pi*lambda_c.*r_dist.^2);
    syms n

    k1_r= (delta/(1-delta)).*Tplot.*r_dist.^2.* hypergeom([1, 1-delta], 2-delta,-Tplot);
    k2_r= symsum((((-1).^(n+1).*((big_X./P_c).^n).*(r_dist.^(2+(n.*alpha))).* hypergeom([1, -delta-n],-delta-n + 1,-Tplot))/factorial(n).*((1+(n*alpha/2)))), n, 0, Inf);
    k3_r= delta.*((P_c./big_X).^delta).*gammainc(delta, (r_dist.^(alpha).*(big_X./P_c)));

    fun = @(r_dist) exp((-lambda_c_eq_intf.*pi).*(k1_r + k2_r -k3_r)).*f_r_dist;
    p_c_uic = integral(fun,0,Inf);
    plot(big_X,p_c_uic)

Answer 1

[Walter Roberson]

Why are you looping over k when you never use k in the loop? You are just doing the same calculation multiple times storing into the same location each time.

Your second integral us over r which would tend to suggest that you should be using r as your variable of integration. But you do not do that: your r is a constant value 5 and your variable of integration is x, which you ignore. Remember that the result of integration of an expression that is constant in the variable of integration is just the expression times the difference in bounds. The difference in bounds is infinite so the result of the integration is going to be sign() of the constant expression times infinity.

When I look at what you are integrating it looks rather likely to me that the constant expression is symbolic and not convertible to numeric as I think I see an unresolved symbol in it. integral() is strictly for numeric results.

[Walter Roberson]

You did not post your complete revised code. I doubt that you have correctly defined the k functions.

[Walter Roberson]

You assign 5 to r and then assign to k1(r). That is an array assignment, not the creation of a function or formula. It defines k1 at exactly one point, when r=5

[Oluwatosin Amodu]

Thank you for those observations and comments. I have taken time to revise the code bearing in mind all you have mentioned. I have removed the code in the comment section and thoroughly revised the one in the question (which is meant to plot only the first equation). I am not still sure what is wrong. Thank you

[Walter Roberson]

You have

      fun = @(r_dist) exp((-lambda_c_eq_intf.*pi).*(k1_r + k2_r -k3_r)).*f_r;

but the expression ignores r_dist . And you are still defining at exactly one r, r = 5.

Look at the image you linked to,

and look at the P_c definition. It involves lambda superscript eq subscript c-intf, which is defined in terms of lambda_c and lambda_d times (P_d/P_c)^delta times nu/mu . But P_c is the equation itself, so P_c is defined recursively unless delta or nu are 0 . Is the model you are dealing with deliberately talking about some kind of equilibrium state?

[Walter Roberson]

The symbolic toolbox is not able to resolve the symsum() in k2_r to numeric values even when vpa() is used, which is why you are getting sym that cannot be converted to numeric for the purposes of numeric integration.

My tests with Maple suggest that the symsum are -infinity for at least some of them. Though I do see some question about whether it should be +infinity or -infinity, with the -infinity seeming more likely.

[Oluwatosin Amodu]

Thank you so much for your follow up. The p_c which equates to the integral is a small letter and the parameter for calculating lambda^eq_c_intf is capital. I also have changed f_r to f_r_dist.

[Walter Roberson]

Your revised

      fun = @(r_dist) exp((-lambda_c_eq_intf.*pi).*(k1_r + k2_r -k3_r)).*f_r_dist;

still ignores r_dist.

Earlier you defined f_r_dist as

      f_r_dist= 2.*pi.*lambda_c.*r_dist.*exp(-pi*lambda_c.*r_dist.^2);

but that uses the value of r_dist as known at the time of the assignment. Which is r_dist = big_X/P_tx_h and those are constant values.

If you want to define f_r_dist in terms of an r_dist to be defined later, then you need to define it as a function, like

f_r_dist = @(r_dist) 2.*pi.*lambda_c.*r_dist.*exp(-pi*lambda_c.*r_dist.^2);

and

        fun = @(r_dist) exp((-lambda_c_eq_intf.*pi).*(k1_r + k2_r -k3_r)).*f_r_dist(r_dist);

[Oluwatosin Amodu]

My objective is to plot the small p_c against big_X. The values of k_r is what I show in the equations. All my details are in the code. I have tried changing to a function as suggested (I have edited the code).

[Walter Roberson]

Is there a P_c different from p_c ? The resolution of the images is not good enough for me tell if there is a p_c or which equations use P_c and which use the hypothetical p_c . It looks to me as if the equations all use P_c and that P_c is defined in terms of itself.

[Oluwatosin Amodu]

Thank you @Walter Roberson. I have made edits to the question for more clarity.

[Jan]

@Oluwatosin Amodu: What does your flag mean? Do you want to ask Walter to remove this image from his comment? If so, please ask him clearly.

[Walter Roberson]

Equations cannot be copyrighted. Images can be copyrighted, but using a single image out of a paper for educational purposes is within fair use / fair dealing exemptions in most countries.