how to implement an infinite series?
5 次查看(过去 30 天)
显示 更早的评论
im trying to implement the function in the attached image but i get that error "
Attempt to reference field of non-structure array.
at the step where (v) is being calculated. i dont know what it means and i dont know how to fix it
if true
clc;
close all;
N=32;
K=14;
c=0.8;
p=1;
pc=p*gamma(1+2/c);
T=399.01923099362968691492508329876;
s=T/(pc*(1+lambda);
syms i j f
nlambda = 50;
lambdavals2 = linspace(0,30,nlambda);
lambdavals = 10.^(lambdavals2/10);
for L = 1 : nlambda
lambda = lambdavals(L);
for i=K:N/2
w=nchoosek(N/2,i) %first series over i
z=arrayfun(@(j)nchoosek(K+i-1,j)*(-1)^(j),0 :K+i-1, 'uniform',0)%second series over j
v=symsum((-1)^(f)*(N+j-K-i+1)^(f)*gamma((f+1)*c/2)/(factorial(f)*(p*(T/(pc*(1+lambda))))^((f+1)*c/2)),i ,0, Inf)%third series over f
pd1temp=w*z*v
end
pd1 = c * K * nchoosek(N/2,K) * pd1temp; %the whole function
end
end
27 个评论
Jan
2017-6-26
编辑:Jan
2017-6-26
No useful comments, outcommented code, no auto-indentation, much code not related to the question, incomplete error message which does not mention, which line is failing: This makes it hard to post an answer. Please invest some time to narrow down the question. Post the relevant part of the code only.
Walter Roberson
2017-6-27
Are there any constraints on the values of c, k, or N ? For example is N constrained to be an even positive integer? Is k constrained to be a non-negative integer less than or equal to N/2 ? Is c constrained to be between 0 and 1?
Fatma Abdullah
2017-6-28
YES, N is constrained to be positive from 0 to 32.....k also constrained to be positive from 0 to 14 ....c could be any positive value from 0 to 2 including fractions
Walter Roberson
2017-6-28
So N could be non-integral or could be odd-valued, so N/2 is not necessarily an integer? Should
(N/2)
( k )
be understood as the generalization
gamma(N/2+1)/(gamma(k+1)*gamma(N/2-k+1))
instead of as
factorial(N/2) / (factorial(k) * factorial(N/2-k))
?
And to confirm, k could be (for example) pi, since 0 < pi < 14?
I notice in your code that you set k = 14. You say that "positive from 0 to 14" and since 0 is not positive, that would suggest that your bounds are "exclusive" rather than "inclusive", the interval (0,14) instead of [0,14] -- but if so then k could not be 14. Please clarify.
Is N required to be an integer? Is N required to be an even integer, which is to say an integer that is a multiple of 2? Is k required to be an integer? Can k be 0? Can c be 0? Can c be 2 ?
Fatma Abdullah
2017-7-2
N is always integer positive even number,so yes N is an integer that is a multiple of 2, k is required to be an integer positive between value>0 and N/2 and could be even or odd so the limits is ]0,N/2] , c is limited between value>0 and infinity including fractions and is required to be positive so the limits of c is ]0,inf[
Walter Roberson
2017-7-4
I have not been able to find any closed form solution for the infinite sum. It looks to me as if numeric approximation might work well, but numeric approximation is not the solution to an infinite summation.
Walter Roberson
2017-7-4
The denominator of your summation includes rho and s to the power of (f+1)*c/2 with f growing to 0. You could imagine that c might be small, but since it is positive, no matter how small it is, (f+1)*c/2 is going to grow indefinitely large as f approaches infinity. So we have a denominator that includes (rho*s) to an indefinitely large power as f approaches infinity.
Now, if rho*s is greater than 1, then as the power (f+1)*c/2 grows indefinitely large, so would (rho*s)^((f+1)*c/2) . That would lead you to an indefinitely large positive denominator. If the numerator of the fraction is growing more slowly then the implication is that the fraction as a whole would converge to 0, allowing for the possibility of a finite summation.
But if rho*s is less than 1, then as the power (f+1)*c/2 grows indefinitely large, then (rho*s)^((f+1)*c/2) would grow indefinitely small. If the numerator of the fraction is not growing smaller at a faster rate, then this would lead to an indefinitely small denominator which would explode, leading to divergence of the sum.
Can (rho*s) be less than 1? If you work through your equations, rho*s = 399.01923099362968691492508329876 / (gamma(2/c + 1)*(lambda + 1)) . Your c is permitted to range up to infinity, so 2/c can decrease indefinitely towards, so gamma(2/c+1) can head towards gamma(1) = 1. So what about lambda? We do not have an expressed limit on lambda so far; looking at the way lambdavals is constructed, it would appear that lambda would always be at least 1 and could be at least 1000, possibly more, so that gives you 399/1000 which is less than 1, so Yes, rho*s can indeed be less than 1. gamma(2/c+1) would have to be less than about 1/4 to prevent that, and that combination cannot happen for non-negative c.
Maybe we had better have another look at the numerator. The N+ etc looks non-negative, so that to the power f is going to grow indefinitely large. gamma((f+1)*c/2) is going to grow indefinitely large. With those two growing indefinitely large, if the denominator grows indefinitely small, we have divergence.
So back to the denominator then, where I have neglected the factorial(f) . Perhaps that could be shown to grow more quickly than (rho*s)^((f+1)*c/2) falls.... ah yes, it appears that for constant rho*s that is indeed the case, that factorial(f) grows faster than (rho*s)^((f+1)*c/2) decreases. So maybe we do not need to worry? But then there is the question of whether the gamma in the numerator cancels that out... I can see circumstances in which it can...
At the moment it is looking to me as if the sum can be divergent, but it needs more examination.
Walter Roberson
2017-7-4
... Yes, with a few test values I can establish that the sum is divergent for c >= 2 and larger lambda > 1, and I can establish that the terms of the sum converge towards 0 for c < 2 and larger lambda. At the moment I am not certain of the behaviour for lambda = 1 exactly.
David Goodmanson
2017-7-5
编辑:David Goodmanson
2017-7-12
Hello fatma and Walter,
I took a look at just the inner, infinite series. For c = 0, the series diverges simply because gamma(0) = inf. Aside from that, and with the help of the large argument expression for the gamma function, for large n the ratio between successive terms is
|A_n+1 / A_n| = B (c/2)^(c/2) n^(c/2 -1)
where
B = (N+j-k-i+1)/(rho s)^(c/2). [modified again, second time]
The ratio test says the series is
divergent for c > 2
convergent for 0 < c < 2
for c = 2:
divergent for |B| > 1
convergent for |B| < 1
undetermined for |B| = 1
Hopefully the numerical results are in accord with this. How quickly the series converges is of course another matter.
David Goodmanson
2017-7-5
编辑:David Goodmanson
2017-7-5
I corrected the comment above to change the expression for B to include ^(c/2), and the convergent-divergent conclusions are the same.
David Goodmanson
2017-7-10
编辑:David Goodmanson
2017-7-10
I wanted to mention that you can replace the infinite sum by a (numerical) integral from 0 to inf. In some cases the sum needs a huge number of terms before the terms start to get smaller, and I think the integral may work over a wider range of parameter values than the sum does.
If you use the identity
gamma(a) = Integral {0,infinity} x^(a-1) e^(-x) dx
on the gamma function in the numerator and sum the series, then add and subtract a known integral, you can arrive at
c2 = c/2;
B = (N+k-j-i+1)/(rho s)^(c/2)
result = (1/(rho*s)^c2)*( (1/(B*c2)) + integral(@(x)b(B,c2,x),0, inf) )
function y = b(B,c2,x)
y = (exp(-x)-1).*x.^(c2-1).*exp(-B*x.^c2);
y(x==0) = 0;
end
qwrqwerth a
Fatma Abdullah
2017-7-11
i think i do understand your contributions; but unfortunately i still haven't been able to understand how to numerically apply this function on matlab
Walter Roberson
2017-7-11
You had indicated that c can be 0 to infinity. We calculate that if c > 2 then the summation is infinite and so does not need to be carried out. We calculate that if c == 2 exactly then under some circumstances the summation is infinite and so does not need to be carried out. Therefor there is no point doing the summation unless c < 2, or unless c == 2 exactly and a particular relationship holds on some of the other variables.
Walter Roberson
2017-7-11
David, you wrote,
B = (N+k-j-i+1)/(rho*s)^(c/2). [modified]
for c = 2: divergent for |B| > 1
That test looks at any one B value, but the B values change because there are two other levels of summation happening that involve changing some of those values. I just wanted to check some things along those lines:
In the (N+k-j-i+1) component, N and k are fixed in all levels of the summation, but i and j can increase relative to their starting values; an increase of i or j would give a smaller numerator, which would lead to a smaller B value (provided the terms are non-negative, which should be the case given the constraints.) Therefore we should be able to test the first combination for B, and if it gives a value < 1 then it is not possible for any of the other possibilities to end up with B > 1.
The starting combination is i = k and j = 0, which would lead to (N+k-0-k+1) = N+1 . So looking overall, if c == 2 and (N+1)/(rho* s)^(c/2) > 1 then there would be divergence. Taking into account that this is the c==2 case, that could be rewritten as N > (rho*s) -1 would give divergence.
Let's see... for c = 2, s works out as T/(rho*(1+lambda)) so rho*s works out as T/(1+lambda), so N > T/(1+lambda) - 1 . The right hand side is least for larger lambda . We are not given a hard constraint on lambda but according to the code it can be up to 10^(30/10) = 1000, and T is about 400. For simplicity let lambda = 999 so 1+lambda = 1000 and approximate T as 400, so the constraint becomes N > 1/4 - 1 so Yes, there would be a problem . The largest supported lambda would become T/(N+1)-1 which would be about 12 for N = 32 and the given T.
On the other hand, the code is evaluating the summation for multiple lambda, so for the given N and T ranges, this mostly tells us that we can stop going fairly early under this circumstance of c = 2 exactly.
Walter Roberson
2017-7-11
The user sent me a copy of the paper, which is "Performance analysis of ordered-statistic greatest of-constant false alarm rate with binary integration for M-sweeps", by X.W. Meng. Glancing at it, two things stand out soon:
- Yes, c=2 is a case of interest; and
- they talk about doing "binary integration", which is not something I am familiar with:
"Binary integration is a non-coherent approach offering most of the benefits of integration with a simpler and less expensive receiver implementation [7]. This technique is based on detecting a specified number (S) of the available pulses (M). If at least S of the pulses are detected, a ‘hit’ is declared. This is sometimes referred to as coincidence or S-out-of-M detection."
I have not read deeply enough into the paper to understand how that affects the situation.
David Goodmanson
2017-7-12
编辑:David Goodmanson
2017-7-12
Hi Walter, thanks for your continued interest in this. Unfortunately I seem to be losing it trying to write the expression for B, which should not be that hard. Per the numerator in the original posting, it's actually
B = (N+j-k-i+1)/(rho s)^(c/2). [modified again]
not
B = (N+k-j-i+1)/(rho s)^(c/2). [modified]
so i and j come in with opposite sign and you can't make the same conclusion that you have. I apologize if this costs you some time.
When c=2 the gamma functions go away and you get a straightforward geometric series that converges for |B|<1. A worse problem is when c < 2 and |B| is only moderately large. This code shows log10 of the individual terms:
n = 0:63000;
c = 1.8;
B = 3;
y = .434*(n*log(B) + gammaln((n+1)*(c/2)) - gammaln(n+1));
plot(n,y)
and it tops out at around n = 23,000 where the value of the term itself is roughly 1e992, eventually dropping to 1e-20 at around n = 62,600. Not a formula for success, so |B| had better be small. I can't vouch for its accuracy but the integral method does give a normal looking answer here, .2847.
Fatma Abdullah
2017-7-12
Walter....binary integration is a method used in the paper and its a definition "binary integration" made to explain a certain method of combining the back scattered signals from a target in the radar receiver , its not a special type of mathematical integration as it might first come to your mind......
Fatma Abdullah
2017-7-12
I've forgot to mention that rho would actually be cancelled as described in the attached img......deeply sorry for not mentioning that earlier
I've tried to code differently as:
if true
clc;
close all;
M=8;
N=32;
K=14;
c=1.2;
pc=p*gamma(1+2/c);
s1=1;
T=54.199575;
syms f
nlambda = 30;
lambdavals2 = linspace(0,30,nlambda);
lambdavals = 10.^(lambdavals2/10);
pd = zeros(1, nlambda);
for L = 1 : nlambda
lambda = lambdavals(L);
Gsum=0;
S=T*p/(p*gamma(1+2/c)*(1+lambda))
pc=p*gamma(1+2/c);
for i=K:N/2
i;
w=nchoosek(N/2,i); %first series over i
z=0;
sumjf1=0;
for j1=0:K+i-1
z=(nchoosek(K+i-1,j1)*((-1)^j1));%second series over j1
v1=subs((((-1)^(f))*((N+j1-K-i+1)^(f))*gamma((f+1)*c/2)/(factorial(f)*((T/(gamma(1+2/c)*(1+lambda)))^((f+1)*c/2)))),f,0,inf);%infinite series over f
v1sum=sum(v1);
pd1temp=v1sum*z
sumjf1=sumjf1+pd1temp;
end
sumjf=sumjf1*w;
Gsum=Gsum+sumjf ;
end
pd1 =1*c* K * nchoosek(N/2,K) * Gsum %the whole function
pdtemp = arrayfun(@(i) nchoosek(M,i)*(pd1)^(i)*(1-pd1)^(M-i), s1:M, 'Uniform' ,0);%pd1 is implemented in this function of pd to draw pd against lambda(db)
pdsum = sum( [pdtemp{:}] );
pd(L) = pdsum;
end
plot(lambdavals2, pd);hold on
axis([0 30 0 1])
end
Fatma Abdullah
2017-7-12
When i draw this equation it was totally different from what i expected( almast a constant straight line with pd=0 for all values of lambda )what i expected is in the attached img (the curve I've pointed)
Walter Roberson
2017-7-12
Your original code had s=T/(pc*(1+lambda); which is missing a ')' and so is ambiguous. The image you just posted has T divided by rho_c all multiplied by 1+lambda but you appear to have implemented as T divided by (rho_c times 1+lambda)
Fatma Abdullah
2017-7-12
编辑:Fatma Abdullah
2017-7-12
its T/[pc*(1+lambda)] and its the same in the image just a different way of writing ....i might have missed the parenthesis in the code
Fatma Abdullah
2017-7-12
and the value of T is not always 399 it could be any value >0 to thousands .....its calculated according to other equations mentioned in the paper
回答(1 个)
Walter Roberson
2017-6-26
You have
i=13;
which assigns a specific numeric value to i . Then, a few lines later, you have
syms pfa1 i T j f
which overwrites i with a symbolic value, equivalent to doing
i = sym('i');
So now i is symbolic with no value other than itself, and has completely lost track of the fact it had previously been assigned 13.
Then later in your code you assume that i is still 13.
7 个评论
Fatma Abdullah
2017-6-26
no i need (i) as a symbolic value actually, i assigned i=13 just to get rid of that warning of considering it complex,its stupid i know, and my code in the question wasn't well written but i've edited it now
Fatma Abdullah
2017-6-26
编辑:Fatma Abdullah
2017-6-26
okay ....i've edited the code to be like that and it takes literally forever,,,is there any way to fix it to run faster and get me results
if true
M=8;
N=32;
K=14;
c=0.8;
SS=1:8;
i=13;
p=1;
pc=p*gamma(1+2/c);
T=399.01923099362968691492508329876
nlambda = 50;
lambdavals2 = linspace(0,30,nlambda);
lambdavals = 10.^(lambdavals2/10);
pd = zeros(1, nlambda);
stop_count = nlambda;
for L = 1 : nlambda
lambda = lambdavals(L);
pd1temp=0;
s=T./pc.*(1+lambda)
for i=K:N/2
w=nchoosek(N/2,i); %first series over i
for j=0:K+i-1
z= nchoosek(K+i-1,j)*(-1)^(j);
% z=arrayfun(@(j)nchoosek(K+i-1,j)*(-1)^(j),0 :K+i-1, 'uniform',0)%second series over j
% v=arrayfun(@(f)
v=symsum((-1)^(f)*(N+j-K-i+1)^(f)*gamma((f+1)*c/2)/(factorial(f)*(p*s)^((f+1)*c/2)),f ,0, Inf);%third series over f
pd1temp1=w*z*v;
pd1temp=pd1temp1+pd1temp;
end
pd1 = c * K * nchoosek(N/2,K) * pd1temp
plot(pd1,lambdavals2);hold on
end
end
end
Torsten
2017-6-27
It does not make sense to use "symsum" if no symbolic variables are defined.
Best wishes
Torsten.
Fatma Abdullah
2017-6-28
it doesnt matter even if i put the line
if true
syms f
end
still nothing changes,,,it takes forever, nothing comes out, and eventually matlab just doesnt respond and shuts down !!!
Walter Roberson
2017-6-29
You asked for an infinite sum; that can potentially take a long time and a lot of memory unless convergence can be established.
That is why I have been asking the questions about the constraints: because knowing the constraints can make it much easier to establish whether something converges or not.
Fatma Abdullah
2017-7-3
N is always integer positive even number,so yes N is an integer that is a multiple of 2, k is required to be an integer positive between value>0 and N/2 and could be even or odd so the limits is ]0,N/2] , c is limited between value>0 and infinity including fractions and is required to be positive so the limits of c is ]0,inf[
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息
标签
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)