数值积分后显示的值范围不正确

Question

我正在尝试使用 MATLAB 绘制以下方程的错误概率，我想使用命令 trapz 进行数值积分，问题是我得到了一个很好的绘图形状，但是值在 y 轴上是错误的，整个曲线应该在 0 和 1.2 之间，但它在 0.492 和 0.5 之间！！谁能告诉我我的代码有什么问题，或者只是给我一个提示？我真的需要帮助。这是我需要绘制的公式（使用 Maketex 编写）：

这是我的代码：

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

您的代码中有一些错误的表达式。另外我怀疑你应该评估循环内的积分。更重要的是，您的整体网格z 似乎太粗糙了。以下代码为我提供了 P(e) 中第二项的 0~1.2 范围

% 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;

我不知道为什么我需要将dz 乘以两次，但这给了我正确的值范围。但我认为这与v 向量值有关。

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

ans =

    0.0002    1.2241