MATLAB: Summing An infinite Series

beginnerfunctionsMATLABseries

Hi !

I'm a little new to MATLAB but am trying to use it to check the values for part of an infinite sum I have to implement in Python.

This is the sum in question.

And this is my code:

function the_sum = p_sum(l,m,n,alpha,beta,gamma)
syms p;
the_sum = syssum(factorial(l+m+n+p+2)/(factorial(l+1+p)*(l+m+2+p)) * (alpha/(alpha + beta + gamma))^p,p,0,Inf);
end

But when I run it I get an error ! I imagine I'm doing something rather simple that is wrong and would be so grateful if someone could point me in the right direction !

Best Answer

I might start by using the correct function. In MATLAB, that is symsum, not syssum. :) These dumb computers, they refuse to read our minds.

You don't tell us what error you got. It is almost always best to show the ENTIRE error message, as that will usually tell what the problem is immediately. Here, it should probably have said something about syssum not existing. But that is only a wild guess on my part, since you also don't show enough to really test what you are doing.

As for doing the computation itself in numerical form, you could probably play around with Stirlings approximation for those factorials. However, symsum seems to handle it well enough.

l = 2;
m = 1;
n = 2;
alp = .3;
bet = .7;
gam = .4;
symsum(factorial(l+m+n+p+2)/(factorial(l+1+p)*(l+m+2+p)) * (alp/(alp + bet + gam))^p,p,0,Inf)
ans =
3757404/14641

Be careful in defining variables named alpha, beta, and gamma, as they are already the names of useful functions in MATLAB. That will cause problems later on.

I would also strongly suggest that you NOT use a variable name of lower case L, thus l. As you can see, it is quite easily mistaken for the number 1, depending on the font used. One day, you will trip over a bug in your code, and not have a clue what is wrong, because you mistyped the variable name as a 1, or the opposite. As such, I STRONGLY recommend you avoid use of the variable names: lower case L (l), upper case i (I), or either case o or O.

Related Solutions

MATLAB: Warning: Unable to find explicit solution for 1 unknown

Not everything you write down has an analytical solution. In fact, most expressions that you can formulate will have no solution at all.

Here alpha is apparently set at 0.9, but n and omega are left as symbolic variables, so completely unknown. Then you used symsum to compute a result, however symsum was unable to give you something terribly useful. So it gave you a bit of a mess, in a form that still had a limit as n goes to infinity. Essentially, MATLAB was unable to eliminate n from the problem. So this really was not strictly a one variable problem in the eyes of MATLAB.

You wish, apparently, to solve for omega, as a function of n. Let me look at what you seem to have done, since you don't show us the code, only the result of the code. So next, I must guess what you did.

alpha = 0.9;
syms n omega
H = symsum(alpha^n *exp(-1i * omega * n),n,0,inf)
H =
piecewise(exp(omega*1i) == 9/10, Inf, exp(omega*1i) ~= 9/10, -(10*exp(omega*1i)*(limit((9/10)^n*exp(-n*omega*1i), n, Inf) - 1))/(10*exp(omega*1i) - 9))

I've temporarily dropped the norm from there in H, as not immediately relevant to anything.

First what does this mean?

This piecewise function in the solution tells us that under some conditions here, there was something "interesting" that happens. That is, when

exp(omega*1i) == 9/10

THEN we get an inf result.

Effectively, it allows us to solve for the case where you would have a problem. Remember that omega can be ANY number at all. ANYTHING. So if we ever have omega as the value

omega = log(9/10)/(1i)

THEN we would get an infinite result. And that sort of makes sense, if I look at the original expression that created this. If the above is true, then we might take the inverse as

exp(-omega*1i) = 10/9

And therefore, if I raise that to any power n, I would see that

exp(-omega*1i*n) = (10/9)^n

So, under that specific condition, for that specific value of omega, the original symsum reduces to the sum of alpha^n*alpha^(-n), from n=0 to inf. but alpha^n*alpha^(-n) is just 1 for all n. So the sum of 1 (a constant) from 0 to infinity must itself diverge to infinity as n gets arbitrarily large.

That is, when omega is specifically log(9/10)/1i, we get complete garbage. But you need to recognize that this condition imples that omega would then be an imaginary number. See that the second half of the piecewise result indicated what happens when omega is ANY other value. Anything at all.

We don't know if omega is a real or complex number. Apparently the symbolic toolbox was not told so. If omega is a real number, then that condition is meaningless, since the problem arises only when omega was a specific imaginary number.

Regardless, the OTHER half of the condition is when omega would have taken on ANY other value. And then, symsum gave up. It effectively told you that the sum was the limit of the sum, as n goes to infinity. As such, MATLAB never actually gave you something useful for the sum. It did not know the answer, if omega was some general value that may be real or complex.

Finally, you ask MATLAB to solve for an ANALYTICAL solution to something it was unable to deal with in the first place. MATLAB cannot solve for the solution of a problem it could not deal with in the first place!

Now, I have no clue what omega might be. Can it be complex? You tell me. But, IF I tell MATLAB that omega was real in the first place, then what happens?

syms n real
syms omega real
H = symsum(alpha^n *exp(-1i * omega * n),n,0,inf)
H =
(10*exp(omega*1i))/(10*exp(omega*1i) - 9)
norm(H)
ans =
(100/abs(10*exp(omega*1i) - 9)^2)^(1/2)

So now something moderately simple drops out. In fact, now I can go back in and put the norm back into the problem. When I wrap a norm around it, it does something, but nothing that exciting. We now have an abs in there now, which does confuse things just a bit. Is that useful? Only you know. I'm just playing around, since you never said what omega means to you.

solve(norm(H) == 0.5)
ans =
Empty sym: 0-by-1
vpasolve(norm(H) == 0.5)
ans =
-1.064710736992428343165281i

So while solve did not wish to solve for omega, such that norm(H) == 0.5, vpasolve was able to do so. At least, vpasolve thinks it has found something. However, what did it find? It found an imaginary result, for what should be the value of omega. And the only was we ever got something out of theis mess was if we effectively excluded complex results in the first place.

That suggests there IS NO REAL SOLUTION to the problem. Maybe we are done? Maybe not.

fplot(norm(H))

And that plot seems to suggest that norm(H) has two local minimizers, that look at least reasonably close to have a value of 0.5, in which case omega would be near 3. So next, we look at the minima of that function.

omfun = matlabFunction(norm(H))
omfun =
  function_handle with value:
    @(omega)sqrt(1.0./abs(exp(omega.*1i).*1.0e+1-9.0).^2.*1.0e+2)
    
[omsol,fval] = fminbnd(omfun,0,5)
omsol =
          3.14159439201346
fval =
         0.526315789473882

And that is also interesting.

fplot(norm(H),[2,4])
grid on
xlabel omega

And that tells us that your problem does in fact appear to have no real solution. The norm(H) never gets quite that small, never actually reaches 0.5. For omega roughly near pi, we found a solution that was close to 0.5, at norm(H)=0.526...

So I guess that leaves MATLAB in the form of an omphaloskeptic, contemplating where to go next. That choice is now up to you.

When you get something strange for a result, look at what it tells you. Think about what you see. There might be useful information content in that answer after all.

MATLAB: Generate beta random number without statistics toolbox

As a followup, I note that the answer by dpb seems to be incorrect.

There are several ways to generate a beta random variable. Simplest is the direct use of betaincinv. (The hardest part of that is spelling betaincinv.) What I do not know is if you can use that utility, in an executable. It should be possible.

For example, let me pick alpha = 2, beta = 3, with a sample size of 1e6.

alpha = 2;
beta = 3;
N = 1e6;
U = rand(N,1);
tic,Z = betaincinv(U,alpha,beta);toc
Elapsed time is 0.432232 seconds.
histogram(Z,100,'norm','pdf');
hold on
H = fplot(@(x) betapdf(x,alpha,beta),[0,1]);
H.Color = 'r';

That is clearly reasonable. So if you can use betaincinv, then you are done. As I show in my comment on the answer by dpb, you can also use the ratio from a couple of gamma random variables. But betaincinv seems easiest. It seems reasonably fast given what it does, at less than .5 seconds for a million such events.

And if you cannot use betaincinv, then gammaincinv would also fail. So as long as you can use betaincinv, use that.

Best Answer

Related Solutions

MATLAB: Warning: Unable to find explicit solution for 1 unknown

MATLAB: Generate beta random number without statistics toolbox

Related Question