How to do numerical integration using Matlab and how to plot it?
조회 수: 4 (최근 30일)
이전 댓글 표시
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:
And here is my Matlab code:
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;
댓글 수: 0
답변 (3개)
John D'Errico
2017년 3월 11일
You cannot use numerical integration on a function that has unknown parameters in it!
Trapz CANNOT be used here. Of course, the fact that your upper limit is infinite is also a problem for a tool like trapz.
Your function has c,v,L,n,q ALL unknown. What do you expect to plot versus what? From your question, I think that even you don't know what you want to plot. Sorry, but there is nothing that you can plot, as long as all of those parameters are unknown.
댓글 수: 16
Walter Roberson
2017년 3월 18일
편집: Walter Roberson
2017년 3월 18일
In your code for C, the portion that is 1 divided by the product of factorials and gammas, you use gamma(M+ks.*.5) in your code, which is a mistake: the formula at that point is for Ks with capital K. This is the only place that Ks with capital K occurs in the formula, so it would have been easy to accidentally use ks (lower-case K) instead.
More to the point: the formula at that point is gamma(M+Ks*Ns/2) and you have missed the Ns
Walter Roberson
2017년 3월 11일
편집: Walter Roberson
2017년 3월 11일
In sufficiently new MATLAB:
syms z
Q = @(v) sym(v, 'r');
%syms c v L n q
c = rand();
v = rand() * c; L = rand() * c; n = rand() * c; q = rand() * c; %"L,n, and q are all indicies of summations specified in c"
Pi = sym('pi');
result = Q(.5) - Q(.5)*c/sqrt(Pi) * vpaintegral(exp(-z*(1+Q(1.5)/v))*z^(L-Q(.5))*1/((z/v+1)^n*(z/v+Q(.5))^q), z, 0, inf)
댓글 수: 10
Walter Roberson
2017년 3월 19일
You are right, I had the wrong sign on exp(-lambdas*Ns/2). The revised equation would be
P__e(nu) = 1/2-(1/2)*(Sum(Sum(Sum(Sum(Sum(Sum(Sum(Sum(Sum((1/4)*nu^(eta+N+1/2)*exp(-(1/2)*lambda__1)*((1/4)*lambda__1^2)^((1/2)*alpha)*((1/4)*lambda__2^2)^((1/2)*beta)*exp(-(1/2)*lambda__2)*((1/4)*lambda__1/(e__m*nu))^eta*((1/4)*lambda__2/(e__m*nu))^N*exp(-(1/2)*lambda__r*N__r)*(1/4)^((1/4)*N__s*k__s-1/2)*exp(-(1/2)*lambda__s*N__s)*(1/4)^((1/4)*N__r*k__r-1/2)*((1/4)*lambda__s*N__s)^M*((1/4)*lambda__r*N__r)^Q*e__m^A*binomial(eta, i)*binomial(N, j)*binomial(N-j, A)*(e__m+1)^(N-j-A)*GAMMA(A+Q+(1/2)*N__r*k__r)*GAMMA(M+i+j+(1/2)*N__s*k__s)*2^(A+Q+(1/2)*N__r*k__r)*(Int(exp(-z*nu*(1+(3/2)/nu))*z^(eta+N-1/2)/((z+1)^(A+Q+(1/2)*N__r*k__r)*(z+1/2)^(M+i+j+(1/2)*N__s*k__s)), z = 0 .. infinity))/(factorial(eta)*factorial(M)*factorial(N)*factorial(Q)*GAMMA(M+(1/2)*N__s*k__s)*GAMMA(Q+(1/2)*N__r*k__r)*GAMMA(beta+N+1)*GAMMA(eta+alpha+1)), A = 0 .. N-j), j = 0 .. N), i = 0 .. eta), Q = 0 .. infinity), M = 0 .. infinity), N = 0 .. infinity), eta = 0 .. infinity), beta = 1-(1/2)*k__2 .. infinity), alpha = 1-(1/2)*k__1 .. infinity))/Pi^(1/2)
with variables
[N__r = 2, N__s = 2, k__1 = 2, k__2 = 2, k__r = 2, k__s = 2, e__m = 1, nu = avg, lambda__1 = 3/10, lambda__2 = 3/10, lambda__r = 1/10, lambda__s = 1/10]
I do not have generated code for this yet.
David Goodmanson
2017년 3월 13일
Hi Deema, (Having looked at and benefited from the comments so far I will take this answer in a slightly different direction). The subject is the integral above, not the multiple for loops that get you there. I'll neglect the (c/2)/sqrt(pi) factor in front.
If you set z = vu, dz = vdu and call the resulting function g, then the integral is
Int f(z) dz
= v^(L+1/2) * Int g(u) du
= v^(L+1/2) * Int exp(-u*(v+3/2)) .* u.^(L-1/2) ./ ((u+1).^n .* (u+1/2).^q) du
and a lot of the v dependence goes out in front. This integral is amenable to numerical integration. How far out you have to integrate depends on the constants but if n and q are >= 0, L <=10 and v is small, if you integrate from u = 0 to 40 for example, by then the integrand is down around e-10 at most and dropping quickly. Larger v and smaller L make things better.
The exception is when L = 0 (if it is an integer) or more generally if ( -1/2 < L < 1/2 ). Then the integrand is infinite at the origin, but still integrable. In that case you can take
h(u) = exp(-u*(v+3/2)).*u.^(L-1/2)
and do the integration as
Int g(u) du = Int (g(u)-h(u)) du + Int h(u) du
The first integral is good at the origin and is done numerically, although you may have to go in and insert the correct g(u)-h(u) = 0 at the origin since Matlab will probably come up with NaN there. The second integral has the analytic solution
Int h(u) du = ((v+3/2)^(-(L+1/2)))*gamma(L+1/2)
and the gamma function is available in Matlab.
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!