How do I do numerical summations or integrations?

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.

Answer 1

[Eric]

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.

Answer 2

[David Goodmanson]

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.

Answer 3

[Natalia Stein]

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

Answer 4

[David Goodmanson]

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.

[Walter Roberson]

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.

Answer 5

[Natalia Stein]

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.

[David Goodmanson]

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.

Answer 6

[Hafizur Rehman]

Solve the model in Hansen(1985) numerically using the following utility function:

U(C,N)=logC+Alog(1-N)

[Rik]

This sounds like a question, not an answer. https://www.mathworks.com/matlabcentral/answers/6200-tutorial-how-to-ask-a-question-on-answers-and-get-a-fast-answer

[Walter Roberson]

How does this relate to the P function being asked about in this Question ?