How do I do numerical summations or integrations?

53 次查看（过去 30 天）
Natalia Stein
Natalia Stein 2017-11-1
评论： Walter Roberson 2021-4-3
I am trying to numerically compute a PDF that includes an infinite summation. In other words it should continue to sum until the summation converges. Here is the equation
where I know s (a vector from 0 to 2 in 0.01 increments) and beta (scalar 0<beta<1). I am using Matlab 2017b, so I do not have MuPad. I also have an integral form of this equation, but it is not a closed form either. I would appreciate any advice on how to set it up.

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

回答（6 个）

Eric
Eric 2017-11-1
The strategy here would be to start at n=1 and keep summing over n until your result no longer changes with each new addition, namely because the new values are insignificantly small, smaller than MATLAB can handle. For example (assuming B==0.3):
p = @(s,B,n) (-1).^(n+1).*gamma(n.*B+1).*sin(n.*B.*pi)./(factorial(n).*s.^(n.*B+1));
n = 1;
s = 0:0.01:2; B = 0.3;
Pold = p(s,B,n);
n = n+1;
Pnew = Pold + p(s,B,n);
while ~isequal(Pold(~isnan(Pold)),Pnew(~isnan(Pnew)))
Pold = Pnew;
n = n+1;
Pnew = Pold + p(s,B,n);
end
The ~isnan condition in the while loop is necessary because situations like p(0,B,1)==Inf and p(0,B,2)=-Inf, so Inf + -Inf = NaN. If you wanted a solution over B as well, set B=(0.01:0.01:0.99)' or something of the sort. If you do not like the NaN(s) in your answer, you should be able to fine-tune the above to your needs. Otherwise, I assume you understand the equation well enough to be able to explain why MATLAB is returning NaN(s) for certain parameter values.
David Goodmanson
David Goodmanson 2017-11-2
Hi Natalia,
What Eric says is true in theory, but as with a lot of these kind of problems,you have to take a look at how many terms are involved. Suppose you want to add terms until they get as small as 1e-6, which seems reasonable. The (-)^n and sin(pi*beta*n) factors are in the range +-1 and can provide some cancellation, but since you are taking no explicit advantage of that and just adding things up as you go along, it's not a bad idea to look at the size of the terms without those factors.
factorial(n) is the same as gamma(n+1). Then for example s = 1, beta = .5
s = 1, beta = .5
n = 1:20;
log10term = log10(gamma(n*beta+1)./(gamma(n+1).*s.^(n*beta+1)));
plot(n,log10term,'o-')
The plot of log10 shows that you only need 13 terms to get below -6, which is not so bad. But things converge just a bit more slowly when beta is close to 1 and s is close to 0.
For large n, the gamma function and the factorial function get large in a hurry so it's necessary to compute their logarithms. Matlab can calculate log(gamma(x)) and using the rule log(s^a) = a*log(s) leads to the expression below.
For the fairly innocuous choice s = .2, beta = .93 then
s = .2, beta = .93
C = 1/log(10); % to convert natural log to log10
n = 1:1e6:1e11;
log10term = C*(gammaln(n*beta+1) - gammaln(n+1) - log(s)*(n*beta+1));
plot(n,log10term)
and the plot shows that you need almost 6e10 terms to get there. That's a lot of terms. If you are in a better range of s and beta you can do all right, but clearly the sum formula has its issues.
Natalia Stein
Natalia Stein 2017-11-2
编辑：Walter Roberson 2017-11-2
Thank you both for the replies. I see now that the summation is really not the way to go! I have the integral version of the function, but that gives me all sort of errors. Here is the function
and here is my set up
p = @(u,t,b) (exp(-(u^b)*cos(pi*b/2)))*(cos((t*u) - (u^b)*sin(pi*b/2)));
b =0.5;
t= [0:0.01:2];
When I try to integrate I get the following message
>> q = integral(p,0,Inf);
Not enough input arguments.
Error in @(u,t,b)(exp(-(u^b)*cos(pi*b/2)))*(cos((t*u)-(u^b)*sin(pi*b/2)))
Error in integralCalc/iterateScalarValued (line 314)
fx = FUN(t);
Error in integralCalc/vadapt (line 132)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);
Error in integralCalc (line 83)
[q,errbnd] = vadapt(@AToInfInvTransform,interval);
Error in integral (line 88)
Q = integralCalc(fun,a,b,opstruct);
>> q = integral(p,0,Inf,'RelTol',1e-8,'AbsTol',1e-13);
Not enough input arguments.
Error in @(u,t,b)(exp(-(u^b)*cos(pi*b/2)))*(cos((t*u)-(u^b)*sin(pi*b/2)))
Error in integralCalc/iterateScalarValued (line 314)
fx = FUN(t);
Error in integralCalc/vadapt (line 132)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);
Error in integralCalc (line 83)
[q,errbnd] = vadapt(@AToInfInvTransform,interval);
Error in integral (line 88)
Q = integralCalc(fun,a,b,opstruct);
I am not sure why.
Should I adjust tolerances?
Is there something wrong with the function?
How can I trouble shoot it?
So far I have assigned discreet values to u, and the function seems to compute them, but why does it not integrate?
Natalia
David Goodmanson
David Goodmanson 2017-11-2
编辑：David Goodmanson 2017-11-2
Hi Natalia,
To use 'integral' you have to make sure that your function can produce a vector output given a vector input of u values. For example u .^ b produces element-by-element powers as required, whereas u ^ b doesn't. Same idea with replacing * with .* . Multiplying by a scalar does not require the extra dot.
b =0.5;
t= 1;
p = @(u,t,b) exp(-(u.^b)*cos(pi*b/2)).*(cos((t*u) - (u.^b)*sin(pi*b/2)));
integral(@(u) p(u,t,b), 0, inf)
I don't think that the integrate function can produce results for a vector of t or b values although it would be interesting if someone weighs in and says that's possible.
  1 个评论
Walter Roberson
Walter Roberson 2017-11-2
For vectors of upper or lower bounds, you have to do something like arrayfun() to handle pairs of bounds one pair at a time.
For integrations with fixed endpoints wanting to produce results for each of a vector or array of coefficients, you can use arrayfun to present the coefficients one-by-one, or you can set 'ArrayValued' to true and process the entire array at once. You generally want to maximize the amount of vectored computation in one go (to within memory limits). I suspect that not using ArrayValued true would be faster for smaller arrays of coefficients, whereas ArrayValued true might be faster when the coefficient arrays are likely to be larger than the number of simultaneous probe points that MATLAB was likely to run through in one batch for the individual version.

请先登录，再进行评论。

Natalia Stein
Natalia Stein 2017-11-2
One more thing! I the integral has the closed form for b=0.5.
but when I try to solve the exponential symbolically using int(p, u), I get the original input back.
  1 个评论
David Goodmanson
David Goodmanson 2017-11-2
编辑：David Goodmanson 2017-11-3
sounds like that is not going to help. The expression can be proven from the original infinite sum though, with the help of a particular identity for the gamma function.
For your last couple of posts (and this one of mine) it's better to post material that is not intended to be an answer as a comment.

请先登录，再进行评论。

Hafizur Rehman
Hafizur Rehman 2021-4-3
Solve the model in Hansen(1985) numerically using the following utility function:
U(C,N)=logC+Alog(1-N)

类别

在 Help Center 和 File Exchange 中查找有关 Matrix Computations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by