Incorrect value range shown after numeric integration

Question

I'm trying to plot the probability of error for the following equation using MATLAB, 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):

Equation

This is my 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;

Seeing how you say the shape is correct but the values are wrong, it seems like you have a problem with normalization (that is, multiplying everything by some wrong constant). You could "cheat" and just change plot(avg,b) to plot(avg,(b-min(b))/(max(b)-min(b))*1.2)... — Dev-iL
@Dev-iL I can't cheat in such a thing!! this is important. i know that the mistake is so tiny, but i really tried everything, i can't find the error. — Layal
Is the integral a part of the sum, or outside the sigma term and uses an independent set of parameters? — Yvon
@Yvon i tried both, but i think it should be inside since it contains eta and N and other variables that are defined inside the for loops. — Layal
You might want to fix the last closing bracket in your TeX code; it shouldn't be superscript and can be omitted entirely even. — Adriaan

Yvon Yvon · Accepted Answer · 2017-03-19T19:01:36

There was a few wrong expressions in your code. Also I suspect you should evaluate the integral within the loop. What is more your integral meshgrid z seems too coarse. The following code gives me a 0~1.2 range for the second term in P(e)

% close all; clear;clc;
Nr=2;Ns=2;
lmda1=.3; lmda2=.3;
lmdas=.1; lmdar=.1;
z= .01:.01:40;
k1=2;k2=2;
kr=2.*Nr;ks=2.*Ns;
ax=0;
avg=z;
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./4).^(alp/2).*(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./((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));
                                C=f1.*(2.^(kr.*.5+Q+A)).*avg.^(eta+N);
                                q2=Q;
                                n2=M+id+jd+ks/2;
                                N2=N;
                                eta2=eta;
                                fun2 = exp(-z.*avg.*(1+1.5./avg)).*z.^(eta2+N2-1./2).*(1./((1+z).^(q2).*(1./2+z).^(n2)));
                                itgrl= trapz(fun2)*.01*.01;
                                v = avg.^(eta2+N2+1./2);
                                ax=ax+v.*C.*itgrl;
                            end
                        end
                    end
                end
            end
        end
    end
end
end
b=.5-.5/pi^.5 *ax;
plot(avg,b);grid;

I don't know why I need to multiply dz twice but this gives me the correct value range. But I think it has something to do with the v vector values.

>> [min(.5/pi^.5 *ax),max(.5/pi^.5 *ax)]

ans =

    0.0002    1.2241

Incorrect value range shown after numeric integration

1 Answers