Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

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

  1. Mar 18, 2017 #1
    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. jcsd
  3. Mar 18, 2017 #2

    jedishrfu

    Staff: Mentor

    I fixed your latex as you needed $$ to bracket the expressions
     
  4. Mar 18, 2017 #3
    great! Thank you.
     
  5. Mar 18, 2017 #4

    jedishrfu

    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.
     
  6. Mar 18, 2017 #5

    lewando

    User Avatar
    Gold Member

    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)
     
  7. Mar 18, 2017 #6
    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: I don't get a reasonable output for my code?
  1. Convergence of my code (Replies: 9)

Loading...