数値積分後に不正確な値の範囲が表示される

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はあまりにも粗すぎるようです。次のコードは、私が二回dzを乗算する必要がある理由私にはわからない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; 

を与えるが、これは私に正しい値の範囲を提供します。しかし、私はそれがvベクトル値と関係があると思います。

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

ans = 

    0.0002 1.2241