MATLAB: How to solve this summation with exponential in Matlab

functionMATLABmatlab functionsum

Best Answer

If you are not careful, you will need to use a tool with literally thusands of digits of precision, as those terms will vary by thousands of digits in magnitude. Worse, if you even try to compute any of those terms, you will find that double precision either underflows or overflows.

vpa(factorial(sym(1247)))
ans =
8.458643820551953682887859707083e+3320

As such, in order to compute that series, you will need to use symbolic computations, or use a tool like my HPF. I'll do it using HPF, although I could have as easily used the symbolc toolbox too.

DefaultNumberOfDigits 100
S = hpf(0);
E = exp(hpf(-1211.5));
T = hpf(1);
for n = 0:1247
  S = S + E*T;
  T = T*1211.5/(n+1);
end
S
S =
0.8494300028975462166085239594512340095256328874564113849900813900015574782646329472378321829350836040

That is the lazy way to solve the problem. It involves no real effort, no thought about how to compute the result.

Can you solve the problem using double precision? Well yes. But again, you MUST be careful. And you must understand how to work with double precision in a way that does not underflow or overflow before you have gotten the result you need. The trick is often to use logs on these problems.

S = 0;
T = -1211.5;
for n = 0:1247
    S = S + exp(T);
    T = T + log(1211.5) - log(n+1);
end
S
S =
         0.849430002892729

In this sum, you should realize that almost every single term is effectively zero. So by using natural logs, I computed the log of every term, then exponentiate only after I have computed the log of the complete term.

As you should see, the two computations agree reasonably well. The double precision term is wrong in the last few digits, but you also need ot consider that every term in that series varies by as much as roughly 3 decimal digits. So that is probably as good as I can easily do using doubles.

Related Solutions

MATLAB: How to create this function on MATLAB

With only 1/k in the denominator, that sum will take forever to give you any degree of precision in the sum. So that degree of precision theta is sort of wasted, because you would need literally trillions of terms to get even double precision accuracy.

Think of it like this: how many terms to do need to get 1000 digit accuracy in a sum with terms that decrease with 1/k?

When k is on the order of 1e16, you are finally approaching double precision, thus eps. But your computer is not fast enough to compute 1e16 terms in any reasonable amount of time. And don't forget that it takes a significant amount of computer time to compute anything in a thousand digits of precision.

But, just imagine that you can compute one term of that sum in a millionth of asecond. (That is being very gracious on my part, too. I just wish I had a computer that was that powerful, but I don't work for the NSA.)

So every second, one million terms. 31 million seconds per yar or so. So 3.1e13 terms per year. Or, roughly 300 years before you compute 1e16 terms. And that only gets you to roughly double precision. But you want to compute all this in 1000 digits. So you will need to compute on the order of 1e1000 terms, before you start to get close to that.

You should be able to do it, just around the time the universe suffers heat death. Of course, by then you will have a better computer.

However, you did ask how you might do this. You can use my HPF toolbox. Or you can use the symbolic toolbox. Done in HPF (you need to download it from the File Exchange first, and install it on your search path.)

DefaultNumberOfDigits 1000
theta = hpf('1.306377883863080690468614492602605712916784585156713644368053759966434053766826598821501403701197395707296960938103086882238861447816353486887133922146194353457871100331881405093575355831932648017213832361522359062218601610856679057215197976095161992952797079925631721527841237130765849112456317518426331056521535131866841550790793723859233522084218420405320517689026025793443008695290636205698968726212274997876664385157661914387728449820775905648255609150041237885247936260880466881540643744253401310736114409413765036437930126767211713103026522838661546668804874760951441079075406984172603473107746775740640078109350834214374426542040853111654904209930908557470583487937577695233363648583054929273872814934167412502732669268404681540626763113223748823800118041206286013841914438857151609189388944789912125543384749359092744422082802260203323027106375022288131064778444817003723336406042118742608383328221769687812353049623008802672211104016065088809718347778314022490821844106377494000232824192701');

Now, write a function to compute the sum. First I need to make sure that I've interpreted that iterated exponent properly.

Do you intend it as theta^(3^n), or or as (theta^3)^n? The difference is of course significant. By default, MATLAB would treat the computation of theta^3^n as the latter form.

Before I go too far, let me see how fast HPF is, working in 1000 digits.

timeit(@() theta^3^n)
ans =
    0.0015156
    
1/ans
ans =
       659.82

So HPF can do 660 such computations per second, at least on my computer.

How fast is the symbolic toolboxhere?

digits 1000
theta = sym('1.306377883863080690468614492602605712916784585156713644368053759966434053766826598821501403701197395707296960938103086882238861447816353486887133922146194353457871100331881405093575355831932648017213832361522359062218601610856679057215197976095161992952797079925631721527841237130765849112456317518426331056521535131866841550790793723859233522084218420405320517689026025793443008695290636205698968726212274997876664385157661914387728449820775905648255609150041237885247936260880466881540643744253401310736114409413765036437930126767211713103026522838661546668804874760951441079075406984172603473107746775740640078109350834214374426542040853111654904209930908557470583487937577695233363648583054929273872814934167412502732669268404681540626763113223748823800118041206286013841914438857151609189388944789912125543384749359092744422082802260203323027106375022288131064778444817003723336406042118742608383328221769687812353049623008802672211104016065088809718347778314022490821844106377494000232824192701');
timeit(@() theta^3^n)
ans =
      0.00138
      
1/ans
ans =
       724.62

So the symbolic TB does 725 such computations per second, thus about the same as HPF. And, yes, my computer is getting a little long in the tooth, but it was quite fast when I bought it some years ago, so it is still decent. Not even close to a mllion per second. I think you are in deep trouble. Really deep. Massively deep. Humongously deep.

Sigh. anyway. I'll use the symbolic toolbox here.

function S = bigsum(n,kterms,theta)
  % precompute this number
  theta3n = theta^3^n;
  
  % initialize the sum. Note that 1/2 is representable exactly by
  % either HPF or by the symbolic toolbox, so I can just subtract 1/2.
  S = theta3n - 1/2;
  
  pie = sym('pi');
  
  for k = 1:kterms
      S = S + sin(2*k*pie*theta3n)/k;
  end
  S = vpa(S/pie);
end

It is not too bad, as long as I restrict the number of terms.

tic,S = bigsum(1,100,theta);toc
Elapsed time is 0.325404 seconds.

So 100 terms in 1/3 of a second. But do we even know the very first digit of that number yet?

for kt = 100:105
vpa(bigsum(1,kt,theta),20)
end
ans =
0.81880906421480201873
ans =
0.82165197521120346388
ans =
0.82334942686820709411
ans =
0.82099366721044812494
ans =
0.8187293340428497687
ans =
0.82046391141319398969

It seems to be converging. SLOWLY. The first digit might be an 8, that is, if this is a convergent series. After a fair amount of time,10000 terms gives:

vpa(bigsum(1,10000,theta),20)
ans =
0.82099823595744050401

Personally, I recommend buying a VERY fast computer. Better yet, you need to reconsider if you can do this computation in a finite amount of time, unless you are willing to wait for a few gazillion years.

MATLAB: Increase the amount of decimals in the hyperbolic functions

You need to understand double precision. So far, I've seen two posts by you that show you are not doing so.

First, your example.

digits 100
vpa(sinh(50))
ans =
2592352764293536022528.0

Now, I'll compute the same result, using my own HPF tool.

sinh(hpf(50,100))
ans =
 2592352764293536232043.726661466742692413734453854430802756139316036287616046268961598412632838901239

As you can see, HPF has no problem. What was the difference? There, I created the number 50, as an HPF number. THEN I computed sinh of that value.

What did you do, and why was it different? You started with the number 50. This is a double precision value. THEN you computed the sinh, as a double precision number, thus limited to about 16 digits of precision. And then, and only then, did you tell MATLAB to convert this to a symbolic result! What should you expect?

Just because you executed the command

digits 100

this does not tell MATLAB that ALL computations will be done in 100 digits of precision. No matter what, MATLAB does double precision computations. digits is a symbolic toolbox command, and only applies to computations done there. But even then, if you pass in a result that is only in double precision, then the symbolic toolbox cannot do magic, and undo what you have done.

So, now try this:

vpa(sinh(sym(50)))
ans =
2592352764293536232043.726661466742692413734453854430802756139316036287616046268961598412632838901239

It looks just like what HPF returned, and for good reason. The difference between the various computations above is important, as well as the reason why HPF and the symbolic TB were able to generate the correct result with no problem.

Finally, be careful with these things. Both HPF and sym had no problem with the number 50, because it is EXACTLY representable as an integer, even in double precision. So you need to understand the difference between the next two statements, and why they have a different result.

vpa(sym(50.12324545))
ans =
50.1232454499999988684066920541226863861083984375
vpa(sym('50.12324545'))
ans =
50.12324545

Best Answer

Related Solutions

MATLAB: How to create this function on MATLAB

MATLAB: Increase the amount of decimals in the hyperbolic functions

Related Question