MATLAB Answers


How can I create this function on MATLAB?

Pablo Arboleya さんによって質問されました 2019 年 3 月 30 日
最新アクティビティ Walter Roberson
さんによって コメントされました 2019 年 3 月 31 日
How can I write this function for the variable "n"? Assuming that I already specified the value of the constant theta.

  18 件のコメント

Pablo Arboleya 2019 年 3 月 30 日
wow, amazing... how can I define this function on MATLAB? I'm not familiar with the language, could you help me?
Walter Roberson
2019 年 3 月 30 日
Using the intepretation that theta^3^n is theta^(3^n) then
theta = sym('1306377883863080690468614492602605712916784585156713644368053759966434053766826598821501403701197395707296960938103086882238861447816353486887133922146194353457871100331881405093575355831932648017213832361522359062218601610856679057215197976095161992952797079925631721527841237130765849112456317518426331056521535131866841550790793723859233522084218420405320517689026025793443008695290636205698968726212274997876664385157661914387728449820775905648255609150041237885247936260880466881540643744253401310736114409413765036437930126767211713103026522838661546668804874760951441079075406984172603473107746775740640078109350834214374426542040853111654904209930908557470583487937577695233363648583054929273872814934167412502732669268404681540626763113223748823800118041206286013841914438857151609189388944789912125543384749359092744422082802260203323027106375022288131064778444817003723336406042118742608383328221769687812353049623008802672211104016065088809718347778314022490821844106377494000232824192701')/sym(10)^1000;
syms n positive
syms pi
result = theta^(3^n) - acot(cot(theta^(3^n)*pi))/pi
Pablo Arboleya 2019 年 3 月 30 日
Thank you!

サインイン to comment.




1 件の回答

John D'Errico
回答者: John D'Errico
2019 年 3 月 30 日
編集済み: John D'Errico
2019 年 3 月 30 日

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 =
ans =
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 =
ans =
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;
S = vpa(S/pie);
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
ans =
ans =
ans =
ans =
ans =
ans =
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:
ans =
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.

  1 件のコメント

Walter Roberson
2019 年 3 月 31 日
>> syms k n theta pi positive
>> -1/2 + theta^(3^n) + symsum(sin(2*k*theta^(3^n)*pi)/k,k,1,inf)/pi
ans =
piecewise(~in(exp(exp(n*log(3))*log(theta)), 'integer'), theta^(3^n) - ((log(1 - exp(-pi*exp(exp(n*log(3))*log(theta))*2i))*1i)/2 - (log(1 - exp(pi*exp(exp(n*log(3))*log(theta))*2i))*1i)/2)/pi - 1/2)
With positive integer n and the known theta value, the condition drops out, leaving a finite expression that with sufficient prodding can be rewritten to something fairly simple.

サインイン to comment.

Translated by