# I don't get a reasonable output for my code?

Tags:
1. Mar 18, 2017

### deema_master

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:

$$P(e)=0.5-\frac{.5}{\sqrt{\pi}}\sum_{\alpha=1-k_1/2}^{\infty}\sum_{\beta=1-k_2/2}^{\infty}\sum_{\eta=0}^{\infty}\sum_{N=0}^{\infty}\sum_{M=0}^{\infty}\sum_{Q=0}^{\infty}\sum_{i=0}^{\eta}\sum_{j=0}^{N}\sum_{A=0}^{N-j}v^{\eta+N+1/2} C \int_0^{\infty}\exp(-z*v(1+1.5/v))z^{\eta+N-.5} \frac{1}{(z+1)^{A+Q+Nrkr/2} (z+.5)^{M+i+j+Nsks/2}} dz$$

where:
$$C=0.25*\exp \left(-\frac{\lambda_1}{2}\right)*\left(\frac{\lambda_1^2}{4}\right)^{\alpha/2}*\left(\frac{\lambda_2^2}{4}\right)^{\beta/2}*\exp \left(-\frac{\lambda_2}{2}\right)*\left(\frac{\lambda_1}{4 *em *v}\right)^{\eta}*\left(\frac{\lambda_2}{4 *em* v}\right)^N*\exp (-\frac{\lambda_r*Nr}{2})*0.25^{Ns*ks/4-0.5}*\exp (-\frac{\lambda_s*Ns}{2})*0.25^{Nr*kr/4-0.5}*(Ns*\lambda_s/4)^M*(Nr*\lambda_r/4)^Q*\frac{1}{\eta! M! N! Q! \Gamma (M+Ks*Ns/2) \Gamma (Q+Nr*kr/2) \Gamma (\beta+N+1) \Gamma (\eta+\alpha+1)}*{em}^A {\eta\choose{i}}{N\choose{j}}{{N-j}\choose{A}}(em+1)^{N-j-A} \Gamma (A+Q+Nr*kr/2) \Gamma (i+j+M+Ns*ks/2)*\left(2^{A+Q+Nr*kr/2)}\right.$$

And here is my Matlab code:
Code (Text):
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;

2. Mar 18, 2017

### Staff: Mentor

I fixed your latex as you needed  to bracket the expressions

3. Mar 18, 2017

### deema_master

great! Thank you.

4. Mar 18, 2017

### Staff: Mentor

This is probably something we can't help you with directly. My suggestion is to pick apart your program and validate that intermediate results are coming out correctly.

Programmers will often break up code into smaller pieces to make it easier to debug and would add print statements to show intermediate results.

5. Mar 18, 2017

### lewando

What is going on with z (and avg)?

You have:
z=0.0001:1:40

Is this what you intend?
z = 0.0001 1.0001 2.0001 3.0001 ...etc... 39.0001 (length of z is 40)

6. Mar 18, 2017

### deema_master

Ye
Yes this is what i want...i know z should go to a larger value...but i get a dimension error so i made their size the same.