4 Apr 2021 11:08:13 UTC

Topic 225153

(moderation:

So, for approx $250 AUD you can get a home license for MatLab including the relevant Australian sales tax @10%. I had a student version about fifteen years ago when I was doing a university computing course, but that lapsed. Now MathWorks has a home/hobbyist category to get their base product. To which one can add, for a modest fee(s), further modules if you take an interest in some direction(s). I think I'll get the Symbolic Math Toolbox soon. Certainly beats paying big bucks for the whole shooting gazoo with all the features, many of which one doesn't care for. Smooth move MathWorks. I have a birthday coming up, and this has been on my bucket list. The other one is Mathematica which has similar licensing plans but roughly twice the price.

Now I've just completed the 'renaissance level' estimate of ** PI** using successive approximations of scribed polygons for the area of a unit circle : 'squeezing' the number out between inner & outer polygon constructs. To be exact some 3*2

One beware : don't run the hardware with too much on it while one is running MatLab. I had BOINC, Spotify, Chrome and Adobe Digital Editions all together when I first pressed the 'run' button for the algorithm with the parameters maxed out. It took a few minutes to get control back via Task Manager ( Windows 7 ) and I think caused one of the shorter E@H work units to time out ie. ** very CPU intensive** even for my eight core machine.

Does anyone else have either MatLab or Mathematica or some other like product ?

Cheers, Mike.

I have made this letter longer than usual because I lack the time to make it shorter ...

... and my other CPU is a Ryzen 5950X :-) Blaise Pascal

Language

Copyright © 2022 Einstein@Home. All rights reserved.

## You might want to check out

)

You might want to check out the free R programming language.

## So, having re-discovered

)

So, having re-discovered numerical integration

, I thought I'd estimate/test the accuracy of MATLAB's intrinsic integrating function^{*}a hard slog divide and conquer algorithm. Here I use f(x) = xvs^{2}and find the area under the curve b/w zero and one :.... analytically the answer is

. You can see how impressively close MATLAB is on that - correct down to 16 decimal places - and more over that intrinsic function accuracy beats a plain vanilla first principles approach that usesone-thirdof the ordinate axis ! Plus MATLAB is way faster than my code. Well I think it's pretty cool anyway .... but I guess that's why MathWorks are in the numerical business. :-)10 billion divisionsCheers, Mike.

* If a function is continuous on a closed interval then one can definitely integrate. Indeed some functions are defined by an integrating process.

I have made this letter longer than usual because I lack the time to make it shorter ...

... and my other CPU is a Ryzen 5950X :-) Blaise Pascal

## OK, this time I integrate

)

OK, this time I integrate based upon the functional values at the centres of ordinate axis divisions. So for sin(x) b/w zero and PI :

.... where the correct analytical sum is

. The intrinsic command is ahead by 4 decimal places on the brute force approach ( even after a billion elements ). I think MATLAB maybe using some variant of the Gaussian quadraturetwoapproach, as it is so fast and way ahead on accuracy.^{*}Cheers, Mike.

One of the top mathematical shortcuts of all time, but a bit of a rabbit warren in detail alas. Gauss : what a really smart guy ... :-)*I have made this letter longer than usual because I lack the time to make it shorter ...

... and my other CPU is a Ryzen 5950X :-) Blaise Pascal

## Now here's a thing :

)

Now here's a thing : integrate

b/w limits ofsqrt(4 - x^{2})andzeroand you have a quarter of a circle of radius two ie. the area is going to betwo. This time I use thePIfor generating sums, which conceptually looks like this :trapezoidal ruleHere the intrinsic function is ahead of the trapezoidal technique by some 5 decimal places, and in turn is some five decimal places from the 'true' value ( out to 20 decimal places ).

Cheers, Mike.

I have made this letter longer than usual because I lack the time to make it shorter ...

... and my other CPU is a Ryzen 5950X :-) Blaise Pascal

## I found this other method

)

I found this other method called

which uses a quadratic/parabolic approximation to the curve for each ordinate division, in detail implementing the so-calledSimpson's Rule:1-4-1 approximation..... now out to 10 billion terms, and that doesn't do much better than the trapezoidal rule with fewer terms for this integrand. Conclusions ? I can come up with these ( one or more ) possibilities:

A - I'm crap at MATLAB coding.

B - I haven't yet found MATLAB's intrinsic technique.

C - the method is knocking on the door of the numeric limits for the MATLAB engine with regard to the basic floating point type.

D - I'm too fussy, after all my code hits the mark to 11 significant figures, and is quite close on the twelfth and thirteenth ( "90" vs "89" ).

Aha! I suspected as much : a small endnote on the Gaussian Quadrature page of Wikipedia says ( referring to a variant of Gauss's technique called

) :Lobatto quadrature"An adaptive variant of this algorithm with 2 interior nodes

^{[3]}is found in GNU Octave and MATLAB ..... "Cheers, Mike.

( Well, with that in hand, I might next investigate the error bounds

iterations curve. See how much pain for how much gain ...... )vsI have made this letter longer than usual because I lack the time to make it shorter ...

... and my other CPU is a Ryzen 5950X :-) Blaise Pascal

## Here's the results of the

)

Here's the results of the same integrand as before - the quarter circle - but do a log/log plot the absolute error

the exponent of the number of ordinate divisions :vs... where the highlighted area indicates that there is no further real gain in accuracy going from a billion to ten billion terms. My guess is that it has hit the limits of the precision of whatever is the default float type that MATLAB uses. So now I know. :-)

Cheers, Mike.

( edit ) Of course, one doesn't quite have to do all that to prove a point. If one enters within the script :

% Good old PI

myPI = 3.14159265358979323846;

and then prints out the value of myPI :

3.141592653589793

11600it has evidently been rounded to the available precision - as witnessed by the highlighted digits - preserving only down to the 15th digit after the decimal point. So 20 places is definitely out of order at least without any other effort. I wonder if MATLAB has variable precision available .... yes it does, in the Symbolic Math Toolbox ... which, perchance, I have just earlier this week purchased ..... the

command. Aha! A shiny new function to play with .....vpaintegralI have made this letter longer than usual because I lack the time to make it shorter ...

... and my other CPU is a Ryzen 5950X :-) Blaise Pascal

## My fascination continues : PI

)

My fascination continues : PI just keeps turning up every where you look, in non obvious ways. Take the function 4/(1 + x

^{2}) and find the area under curve b/w zero and one ( ie. integrate with those limits ) then :Not evidently/geometrically PI-ish, but PI sneaks in because the anti-derivative of the above function is the inverse tangent function, where tan

^{-1}(1) = PI/4 ....Cheers, Mike.

I have made this letter longer than usual because I lack the time to make it shorter ...

... and my other CPU is a Ryzen 5950X :-) Blaise Pascal

## The function 1/(1 + x2) of

)

The function 1/(1 + x

^{2}) of itself has some fascination too, as it can be represented by the infinite power series { expand using the binomial theorem or just long divide } :1/(1 + x

^{2}) = 1 - x^{2}+ x^{4}- x^{6}+ x^{8}+ .... + (-1)^{k}x^{2k }+ ......... I say represented, but the above equality only truly holds if the right hand side converges to a limit, which it does if x < 1 { if x = 1 then it oscillates b/w zero and one forever, and diverges if x > 1 }.

Now if one integrates the right hand side term by term { it's just an infinite polynomial so each x

^{2k}integrates to x^{2k+1}/(2k+1) } giving :Integral[(1 + x

^{2})^{-1}] = x - x^{3}/3 + x^{5}/5 - x^{7}/7 + ........ logically this ought be the series representation of tan

^{-1}(x). Which it is for x < 1. What of x = 1 ?1 - 1/3 + 1/5 - 1/7 + 1/9 .......

.... while zig-zag over-shooting and under-shooting this does converge to tan

^{-1}(1) = pi/4 as expected, though really slowly eg. hundreds of terms to get a couple of decimal places of accuracy. Still doesn'twhy the ratio of a circle's circumference to it's diameter is involved though ! What geometric construct, involving circles, yields (1 + xexplain^{2})^{-1}as a solution ? Roger Penrose in his epic book{ chapter fourThe Road To RealityMagical Complex Numbers} has an interesting take on all this, as he examines (1 + z^{2})^{-1}in the complex plane rather than merely (1 + x^{2})^{-1}on the real line.Cheers, Mike.

( edit ) OK. Try this :

.... an actual tangent line to the circle whose radius is one, where (1 + x

^{2})^{-1}is the square of the cosine of the angle whose tangent is x .... :-)I have made this letter longer than usual because I lack the time to make it shorter ...

... and my other CPU is a Ryzen 5950X :-) Blaise Pascal

## The infinite series with the

)

The infinite series with the partial sums :

S

_{n}= 1 + 1/2 + 1/3 + 1/4 + ..... + 1/ndoes not converge, as it is unbounded and can exceed as big a sum as you like just by including enough terms. However the infinite series with the partial sums :

S

_{n}= 1 + 1/4 + 1/9 + 1/16 + .... + 1/n^{2}does converge to the value of PI

^{2}/6 ( to be more precise, you can get as close as you like to that number by adding in more terms, but without exceeding that value ). The best explanation I have seen for this rather obtuse result is on's YouTube channel 3Blue1Brown. In the presentation called Why is pi here? And why is it squared? A geometric answer to the Basel problem he gives a reasonably lucid account of how it comes about. He uses the inverse square law of light as radiation, of all things, plus the Inverse Pythagorean Theorem ( yes, that was news to me too ). But dare I say that the originators of the problem didn't think along such lines. My MATLAB code gives :Grant Sandersonas it's best estimate on my machine before something breaks. Proving it in detail requires the use of

.complex numbersOne could reasonably ask the question : for what value of

doespS

_{n}= 1 + 1/2^{p}+ 1/3^{p}+ 1/4^{p}+ ..... + 1/n^{p}such a series of partial sums converge ? Due to a method called the Comparison Test the series with 1/n

^{3}as terms converges whereas 1/n^{0.5}doesn't. As above it doesn't for= 1 and it does forp=2, so maybe some value of p in between is the limiting case ? It turns out that for anyp> 1 ( ie. strictly greater than one ) the series will converge, however slowly. This requires some high powered concepts like the Maclaurin Integral Test to prove though.pNext up there is the product by

which converges to PI/2 as more multiplicative factors are included, and invented long before either Calculus or the Binomial Theorem were published :Mr JohnWallisPI/2 = (2/1)*(2/3)*(4/3)*(4/5)*(6/5)*(6/7)* ..... *(2n/2n-1)*(2n/2n+1) .....

which was serious propeller head stuff then, and still is now, and even more amazing to think how he went about proving that back in 1656. My MATLAB gives :

Mr Wallis even has a sequence that converges to the square root of PI, also very-bloody-non-obviously as the limit of the following formula when the natural number

tends to infinity :n((n!)

^{2}*(2^{2n}))/((2n)!*sqrt(n))Where the

character indicates that the factorial function is being used { n factorial = the product of all the integers from n down to 1 = n*(n-1)*(n-2)*.... *3*2*1 }. Arrgh ! Again via MATLAB :!There were some seriously clever and dedicated people about. It was detail like those above that led to modern mathematics and it's 'unreasonable effectiveness' in describing reality via physics and other such sciences. These were renaissance answers to classical questions like : what method, by using only a ruler and a compass, yields a square of area the same as a given circle ? This requires, say, knowledge of the numerical value of the square root of PI, which could thus be the side length of that square if the circle was of radius 1.

Cheers, Mike.

( edit ) 'Ruler and compass' is a coded form of : try it and see if you can by a rigorous geometric method. The point is that in this case you can't, as PI is not simply an irrational number but a transcendental one. Rational numbers are ratios of integer lengths, hence the name. Irrational numbers aren't. In a certain technical sense 'most' irrational numbers are transcendental

^{*}: that means they not the solution to any polynomial equation with integer coefficients..... and by extension that in turn means most real numbers are transcendental. So they are an important class of numbers & includes other luminaries such as the natural logarithm base (*) and various combinations likee- yes, that ise^{PI}to the power ofe. But enough said for tonite methinks. :-)PII have made this letter longer than usual because I lack the time to make it shorter ...

... and my other CPU is a Ryzen 5950X :-) Blaise Pascal

## Now here's thing. Look at

)

Now here's a thing. Look at these partial sums again :

S

_{n}= 1 + 1/2^{p}+ 1/3^{p}+ 1/4^{p}+ ..... + 1/n^{p}and now study in more detail the dependence of the final sum upon

, forpvalues really, reallypto oneclose. By the previous post these should always give convergent series. You get :but just above oneWhere the data are plotted for successive runs of the summating machine. The first point on the bottom left is the PI

^{2}/6 limit. Then are plotted the ultimate sums after 100000000 terms each, for values ofcloser and closer to 1 eg. the last point on the top right is for p = 1 + 10p^{-15}. So the pattern revealed is a limit of limits as it were. A meta-limit. Now anyone can see the likely limiting value is. An interesting quasi-fact. This however is proof of nothing at all, as I'm just experimenting. I could be merely disclosing some regularity of the MATLAB code or implementation. But I could still ask what is19doing here ? It's a prime number ....19What I'm getting at is that one can

numbers by limiting processes.defineCheers, Mike.

I have made this letter longer than usual because I lack the time to make it shorter ...

... and my other CPU is a Ryzen 5950X :-) Blaise Pascal