數值積分後顯示的數值範圍不正確

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

有代碼中的一些錯誤的表達。另外我懷疑你應該評估循環中的積分。什麼是你的整體meshgrid z似乎太粗糙。下面的代碼給了我一個0~1.2範圍爲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; 

我不知道爲什麼我要乘dz兩次，但這個給了我正確的值範圍。但我認爲它與v矢量值有關。

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

ans = 

    0.0002 1.2241