How to do numerical integration using Matlab and how to plot it?

I'm trying to plot the probability of error for the following equation using Matlab software, i want to use the command "trapz" for the numerical integration, the problem is that i get a fine shape for the plot, but the values in the y axis are wrong, the whole curve should be between 0 and 1.2 but it is between 0.492 and 0.5!! Can anyone just tell me what is wrong in my code, or just give me a hint? I really need help. Here is my formula that i need to plot written using Maketex:

And here is my Matlab code:

    close all; clear;clc;
    Nr=2;Ns=2;
    lmda1=.3; lmda2=.3;
    lmdas=.1; lmdar=.1;

    z= 0.0001:1:40;
    k1=2;k2=2;
    kr=2.*Nr;ks=2.*Ns;
    ax=0;
    avg=0.0001:1:40;
    em=1;
    ch=2;

    for alp=1-k1.*.5:ch
    for beta=1-k2.*.5:ch
        for eta=0:ch
            for N=0:ch
                for M=0:ch
                    for Q=0:ch
                        for id=0:eta
                            for jd=0:N
                                for A=0:N-jd
                                    %

                                    up=.25.*exp(-lmda1./2).*(lmda1./2).^(alp).*(lmda2.^2./(4)).^(beta./2).*exp(-lmda2./2).*(lmda1./(4.*em.*avg)).^eta.*(lmda2./(4.*em.*avg)).^N.*exp(-lmdas.*Ns.*.5).*.25.^(ks.*.25-.5).*exp(-lmdar.*Nr.*.5).*.25.^(kr.*.25-.5).*(Ns.*lmdas.*.25).^M.*(Nr.*lmdar.*.25).^Q;

                                    cy=up.*(1./(factorial(eta).*factorial(N).*factorial(M).*factorial(Q).*gamma(eta+alp+1).*gamma(N+beta+1).*gamma(M+ks.*.5).*gamma(Q+kr.*.5)));
                                    cj=cy.*(factorial(eta)./(factorial(id).*factorial(eta-id))).*(factorial(N)./(factorial(jd).*factorial(N-jd))).*gamma(M+id+jd+ks.*.5);
                                    f1=(cj.*(factorial(N-jd)./(factorial(A).*factorial(N-jd-A))).*em.^A.*(((em+1).^(N-jd-A))).*gamma(kr.*.5+Q+A));

                                    f2=f1.*(2.^(kr.*.5+Q+A)).*avg.^(eta+N);
                                    ax=ax+f2;

                                end
                            end
                        end
                    end
                end
            end
        end
    end
end
q2=2;n2=2;N2=1;eta2=1;

    fun2 = exp(-z.*avg.*(1+1.5./avg)).*z.^(eta2+N2-1./2).*(1./((1+z).^(q2).*(1./2+z).^(n2)));

     out= trapz(z,fun2);

    b=.5.*(1-ax.*(1./sqrt(pi)).*out.*avg.^(1./2));
    plot(avg,b);grid;

Answer 1

[John D'Errico]

You cannot use numerical integration on a function that has unknown parameters in it!

Trapz CANNOT be used here. Of course, the fact that your upper limit is infinite is also a problem for a tool like trapz.

Your function has c,v,L,n,q ALL unknown. What do you expect to plot versus what? From your question, I think that even you don't know what you want to plot. Sorry, but there is nothing that you can plot, as long as all of those parameters are unknown.

[Deema42]

L,n, and q are all indicies of summations specified in c, so they are known! v can take any positive value. by plotting i mean, i need to plot the result vs something, before i think of doing numerical integration, i tried solving the integration, and the result should be plotted versus v in the equation, since for sure the variable z will not appear anymore, but since it's the first time i do numerical integration, i don't know the result should be versus what in the plot! instead of trapz, what can i use to do numerical integration?

[Roger Stafford]

That first sentence, “L,n, and q are all indicies (sic) of summations specified in c, so they are known!”, would seem to hint that L, n, and q are integers and that some kind of summation is involved. Before you can receive any useful kind of help, you need to greatly expand on that particular sentence.

1) In what way are L, n, and q specified in c - details?

2) Is the summation to occur before integration or afterward - in other words, are there to be many integrations or just one?

3) Can you carefully restate your problem in such a way that your desired result depends only on the parameter, v?

Well-thought-out answers to all these questions (and perhaps a few more) will probably enable either you or one of us to figure out what final result you wish to obtain and what it is to be plotted against.

[Deema42]

Here is my code and the use of trapz is also included :

    close all; clear;clc;
    Nr=2;Ns=2;
     lmda1=.3; lmda2=.3;
     lmdas=.1; lmdar=.1;v=1;
     z= 0.1:.1:10000;
     k1=2;k2=2;
       kr=2.*Nr;ks=2.*Ns;

     ax=0;

     avgdb=0:1:40;
     avg=10.^(avgdb./10);
     em=1;
      ch=2;
      for alp=1-k1.*.5:ch
    for beta=1-k2.*.5:ch
        for eta=0:ch
            for N=0:ch
                for M=0:ch
                    for Q=0:ch
                        for id=0:eta
                            for jd=0:N
                                for A=0:N-jd

                                    qx=A+Q+kr.*.5;
                                    n=M+id+jd+ks.*.5;

                                    up=.25.*exp(-lmda1./2).*(lmda1./2).^(alp).*(lmda2.^2./(4)).^(beta./2).*exp(-lmda2./2).*(lmda1./(4.*em.*avg)).^eta.*(lmda2./(4.*em.*avg)).^N.*exp(-lmdas.*Ns.*.5).*.25.^(ks.*.25-.5).*exp(-lmdar.*Nr.*.5).*.25.^(kr.*.25-.5).*(Ns.*lmdas.*.25).^M.*(Nr.*lmdar.*.25).^Q;

                                    cy=up.*(1./(factorial(eta).*factorial(N).*factorial(M).*factorial(Q).*gamma(eta+alp+1).*gamma(N+beta+1).*gamma(M+ks.*.5).*gamma(Q+kr.*.5)));
                                    cj=cy.*(factorial(eta)./(factorial(id).*factorial(eta-id))).*(factorial(N)./(factorial(jd).*factorial(N-jd))).*gamma(M+id+jd+ks.*.5);
                                    f1=(cj.*(factorial(N-jd)./(factorial(A).*factorial(N-jd-A))).*em.^A.*(((em+1).^(N-jd-A))).*gamma(kr.*.5+Q+A));

                                    %

                                    fun2 = exp(-z.*(1+1.5./avg)).*z.^(eta+N-1./2).*(1./((1+z./avg).^(qx).*(1./2+z./avg).^(n)));

                                    out= trapz(z,fun2);

                                    f2=f1.*(1./2.^(kr.*.5+Q+A)).*out;
                                    ax=ax+f2;

                                    %
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

    b=.5.*(1-ax.*(1./sqrt(pi)));
     semilogy(avgdb,b);

c is these 9 summations, avg is the v in the integral image in my question, but it should be what i want to plot the function with, so i know you will get an error of dimension if you run my code, what should i put this avg in the equation related to the integral, which is "fun2", what should i do in a case like this?

[Walter Roberson]

In your line

fun2 = exp(-z.*(1+1.5./avg)).*z.^(eta+N-1./2).*(1./((1+z./avg).^(qx).*(1./2+z./avg).^(n)));

your z is length 100000 and your avg is length 41. What size of result are you expecting?

Perhaps you should be using ndgrid() to construct matrices of z and avg to operate over, or perhaps you should change the orientation of one of the two and use bsxfun(), or perhaps you should change the orientation of one of the two and use the enhancements from R2016b onward that allow explicit bxsfun() to be omitted.

[Deema42]

since i want to plot "b" versus avg, so i want the size of the result to be the same as avg, z needs to be very large since it should be going towards infinity, what exactly should i do?

[Deema42]

i want to ask something, is it right that the result of trapz "out" is inside the loops as i made in my code??

[Walter Roberson]

You trapz against z. Do you perhaps mean that you want out to be the same size as avg?

[Deema42]

out can't be the same as avg in size, since z is going to infinity, right?

[Walter Roberson]

I have no idea how the code you posted relates to the formula you posted. The only variable they have in common is z, and maybe Q is the same as q. But your code appears to be doing a 9-dimensional integration plus another dimension for avg that is supposed to fit in some-how. We are lost.

[Deema42]

I wrote other letters for L and q just for abbreviation, but it's the same formula...c in the image is a series of 9 summations (9 for loops) and it contains the constants i have in f1...then the integral comes wich is fun2... And in b is the total formula.

[Walter Roberson]

In your image, c acts as a constant multiplier on the integral, so the integral should be outside of the 9 loops.

[Deema42]

Okay but the integral contains (eta,N,n,qx) that are defined in the for loops...so it should be inside it right??

[Walter Roberson]

Ah, you have changed your formula rather a lot!

[Deema42]

yes i wrote the whole formula, the shape of my final plot is fine but the values in the y axis are not reasonable, i don't know where i made a mistake!

[Walter Roberson]

In your code for C, the portion that is 1 divided by the product of factorials and gammas, you use gamma(M+ks.*.5) in your code, which is a mistake: the formula at that point is for Ks with capital K. This is the only place that Ks with capital K occurs in the formula, so it would have been easy to accidentally use ks (lower-case K) instead.

More to the point: the formula at that point is gamma(M+Ks*Ns/2) and you have missed the Ns

[Deema42]

Yes it's lower case k. And i defined ks and kr at the beginning of the code as they are multiplied by Ns and Nr respectively, so it's the same.

Answer 2

[Walter Roberson]

In sufficiently new MATLAB:

syms z
Q = @(v) sym(v, 'r');

%syms c v L n q
c = rand();
v = rand() * c; L = rand() * c; n = rand() * c; q = rand() * c;   %"L,n, and q are all indicies of summations specified in c"

Pi = sym('pi');
result = Q(.5) - Q(.5)*c/sqrt(Pi) * vpaintegral(exp(-z*(1+Q(1.5)/v))*z^(L-Q(.5))*1/((z/v+1)^n*(z/v+Q(.5))^q), z, 0, inf)

[Deema42]

what is vpaint exactly? i don't have it on my Matlab, the version is 2014

[Walter Roberson]

Sorry, that was a typo, it should have been vpaintegral(). Also, I accidentally left out my definition for Q; I have amended the post.

https://www.mathworks.com/help/symbolic/vpaintegral.html

However, that requires R2016b or later.

The closest equivalent before that was

syms z
Q = @(v) sym(v, 'r');

%syms c v L n q
c = rand();
v = rand() * c; L = rand() * c; n = rand() * c; q = rand() * c;   %"L,n, and q are all indicies of summations specified in c"

Pi = sym('pi');
result_sym = Q(.5) - Q(.5)*c/sqrt(Pi) * int(exp(-z*(1+Q(1.5)/v))*z^(L-Q(.5))*1/((z/v+1)^n*(z/v+Q(.5))^q), z, 0, inf);
 result = vpa(result_sym);

However, it appears that the integral is not sufficiently well behaved for the numeric integration from the symbolic package; the result comes back in terms of an internal MuPAD numeric::int call.

If you are willing to integrate starting from 1E-24 instead of 0 then vpa() with the default digits setting can handle that. If you start from 1E-25 then vpa(result_sym,50) can handle that. Below that, such as 1E-26, is a struggle.

This in turn implies that you will probably have problems if you use numeric integration with this integral.

The integral has a very step rise near 0 -- a rise that gets steeper near 0. A quick estimate is that near 10^(-N), the slope approaches 10^(0.97*N) -- so for example near 10^(-10000) the slope is 2.019338199*10^9713 for the random values I happened to substitute. This is not something a numeric integral is going to be happy with.

[Walter Roberson]

Please double-check my entering of the formula:

P(e) = 1/2-(1/2)*(Sum(Sum(Sum(Sum(Sum(Sum(Sum(Sum(Sum((1/4)*nu^(eta+N+1/2)*exp(-(1/2)*lambda__1)*((1/4)*lambda__1^2)^((1/2)*alpha)*((1/4)*lambda__2^2)^((1/2)*beta)*exp(-(1/2)*lambda__2)*((1/4)*lambda__1/(e__m*nu))^eta*((1/4)*lambda__2/(e__m*nu))^N*exp(-(1/2)*lambda__r*N__r)*(1/4)^((1/4)*N__s*k__s-1/2)*exp((1/2)*lambda__s*N__s)*(1/4)^((1/4)*N__r*k__r-1/2)*((1/4)*lambda__s*N__s)^M*((1/4)*lambda__r*N__r)^Q*e__m^A*binomial(eta, i)*binomial(N, j)*binomial(N-j, A)*(e__m+1)^(N-j-A)*GAMMA(A+Q+(1/2)*N__r*k__r)*GAMMA(M+i+j+(1/2)*N__s*k__s)*2^(A+Q+(1/2)*N__r*k__r)*(int(exp(-z*nu*(1+(3/2)/nu))*z^(eta+N-1/2)/((z+1)^(A+Q+(1/2)*N__r*k__r)*(z+1/2)^(M+i+j+(1/2)*N__s*k__s)), z = 0 .. infinity))/(factorial(eta)*factorial(M)*factorial(N)*factorial(Q)*GAMMA(M+(1/2)*N__s*k__s)*GAMMA(Q+(1/2)*N__r*k__r)*GAMMA(beta+N+1)*GAMMA(eta+alpha+1)), A = 0 .. N-j), j = 0 .. N), i = 0 .. eta), Q = 0 .. infinity), M = 0 .. infinity), N = 0 .. infinity), eta = 0 .. infinity), beta = 1-(1/2)*k__2 .. infinity), alpha = 1-(1/2)*k__1 .. infinity))/Pi^(1/2)

[Walter Roberson]

From your code, I deduce that:

[N__r = 2, N__s = 2, k__1 = 2, k__2 = 2, k__r = 4, k__s = 4, e__m = 1, nu = avg, lambda__1 = 3/10, lambda__2 = 3/10, lambda__r = 1/10, lambda__s = 1/10]

where you define

avg = 0.0001:1:40

at one point. That vector definition does not fit in with the rest of the equation.

I also notice that your left hand side, P(e) involves e, which otherwise does not appear.

You speak of a y axis, but what is that y axis, and what is your x axis? Your only non-scalar variable that is not bound in by Sum() or Int() is avg . Could it be that your definition should be for P(avg) instead of P(e) -- that is, that your x axis is the avg vector in the range 0 to 40, and that your y axis is the resulting P value?

[Deema42]

It'a an abbreviation for probability of error...yes the x axis is the avg, and y axis is the result of the function evaluated at avg.

[Deema42]

Why the variable avg doesn't fit?

[Deema42]

One of the exponentials should be powered by (-lambdas*Ns/2) in your code...

[Walter Roberson]

avg did not fit because it was a vector when everything else was a scalar and it appeared that the result was supposed to be a scalar. With the information that avg is the independent variable then things work out for the formula.

I will recheck the lambdas later in the day.

[Deema42]

Okay...how can i fix this?

[Walter Roberson]

You are right, I had the wrong sign on exp(-lambdas*Ns/2). The revised equation would be

P__e(nu) = 1/2-(1/2)*(Sum(Sum(Sum(Sum(Sum(Sum(Sum(Sum(Sum((1/4)*nu^(eta+N+1/2)*exp(-(1/2)*lambda__1)*((1/4)*lambda__1^2)^((1/2)*alpha)*((1/4)*lambda__2^2)^((1/2)*beta)*exp(-(1/2)*lambda__2)*((1/4)*lambda__1/(e__m*nu))^eta*((1/4)*lambda__2/(e__m*nu))^N*exp(-(1/2)*lambda__r*N__r)*(1/4)^((1/4)*N__s*k__s-1/2)*exp(-(1/2)*lambda__s*N__s)*(1/4)^((1/4)*N__r*k__r-1/2)*((1/4)*lambda__s*N__s)^M*((1/4)*lambda__r*N__r)^Q*e__m^A*binomial(eta, i)*binomial(N, j)*binomial(N-j, A)*(e__m+1)^(N-j-A)*GAMMA(A+Q+(1/2)*N__r*k__r)*GAMMA(M+i+j+(1/2)*N__s*k__s)*2^(A+Q+(1/2)*N__r*k__r)*(Int(exp(-z*nu*(1+(3/2)/nu))*z^(eta+N-1/2)/((z+1)^(A+Q+(1/2)*N__r*k__r)*(z+1/2)^(M+i+j+(1/2)*N__s*k__s)), z = 0 .. infinity))/(factorial(eta)*factorial(M)*factorial(N)*factorial(Q)*GAMMA(M+(1/2)*N__s*k__s)*GAMMA(Q+(1/2)*N__r*k__r)*GAMMA(beta+N+1)*GAMMA(eta+alpha+1)), A = 0 .. N-j), j = 0 .. N), i = 0 .. eta), Q = 0 .. infinity), M = 0 .. infinity), N = 0 .. infinity), eta = 0 .. infinity), beta = 1-(1/2)*k__2 .. infinity), alpha = 1-(1/2)*k__1 .. infinity))/Pi^(1/2)

with variables

[N__r = 2, N__s = 2, k__1 = 2, k__2 = 2, k__r = 2, k__s = 2, e__m = 1, nu = avg, lambda__1 = 3/10, lambda__2 = 3/10, lambda__r = 1/10, lambda__s = 1/10]

I do not have generated code for this yet.

Answer 3

[David Goodmanson]

Hi Deema, (Having looked at and benefited from the comments so far I will take this answer in a slightly different direction). The subject is the integral above, not the multiple for loops that get you there. I'll neglect the (c/2)/sqrt(pi) factor in front.

If you set z = vu, dz = vdu and call the resulting function g, then the integral is

Int f(z) dz
=  v^(L+1/2) * Int g(u) du
=  v^(L+1/2) * Int exp(-u*(v+3/2)) .* u.^(L-1/2) ./ ((u+1).^n .* (u+1/2).^q) du

and a lot of the v dependence goes out in front. This integral is amenable to numerical integration. How far out you have to integrate depends on the constants but if n and q are >= 0, L <=10 and v is small, if you integrate from u = 0 to 40 for example, by then the integrand is down around e-10 at most and dropping quickly. Larger v and smaller L make things better.

The exception is when L = 0 (if it is an integer) or more generally if ( -1/2 < L < 1/2 ). Then the integrand is infinite at the origin, but still integrable. In that case you can take

    h(u) = exp(-u*(v+3/2)).*u.^(L-1/2)

and do the integration as

Int g(u) du   =   Int (g(u)-h(u)) du  + Int h(u) du

The first integral is good at the origin and is done numerically, although you may have to go in and insert the correct g(u)-h(u) = 0 at the origin since Matlab will probably come up with NaN there. The second integral has the analytic solution

Int h(u) du = ((v+3/2)^(-(L+1/2)))*gamma(L+1/2)

and the gamma function is available in Matlab.

[Deema42]

it is a good idea to make the variable v outside by making uv=z, the shape is fine for the error probability but the values on the y axis are not reasonable, they should be between 0 and 1.2 , but thet are between .48 and .5 which is wrong, can you suggest a solution to this or at least to give me a hint, i made an update in my question. Thank you.