How to do numerical integration using Matlab and how to plot it?

2 次查看（过去 30 天）
Deema42
Deema42 2017-3-11
评论： Walter Roberson 2017-3-19
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;

请先登录，再回答此问题。

回答（3 个）

John D'Errico
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 个评论
显示 14更早的评论
Walter Roberson
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
Deema42
Deema42 2017-3-19
Yes it's lower case k. And i defined ks and kr at the beginning of the code as they are multiplied by Ns and Nr respectively, so it's the same.

请先登录，再进行评论。

Walter Roberson
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 个评论
显示 8更早的评论
Walter Roberson
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
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.
  1 个评论
Deema42
Deema42 2017-3-18
it is a good idea to make the variable v outside by making uv=z, the shape is fine for the error probability but the values on the y axis are not reasonable, they should be between 0 and 1.2 , but thet are between .48 and .5 which is wrong, can you suggest a solution to this or at least to give me a hint, i made an update in my question. Thank you.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by