I'm A Shameless Exponent of MatLab (LOL!)
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*21000 polygons and PI correct to the book value of some 20 decimal places using The Euclidean Algorithm ( which is independent of any knowledge of PI in advance ). I know, a mere 20 places, but so quickly ! Once upon a time there was this Dutch guy who spent half a lifetime doing it. I'll attack it with Newton's approach next, he blew the doors off previous attempts. I've also got a good estimate of e - the base for natural logarithms - but with the much slower converging (1+1/n)n sequence. How strange, in a way, to have numbers only formally defined with convergence sequences. Back in their day computers were people who did, well, computations. I stand in awe of say, Kepler who did all his great works before the invention of logarithms. I'm using A First Course In Mathematical Analysis by David Alexander Brannan as a theory guide.
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
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 vs a hard slog divide and conquer algorithm. Here I use f(x) = x2 and find the area under the curve b/w zero and one :
.... analytically the answer is one-third. 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 uses 10 billion divisions of 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. :-)
Cheers, 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 two. 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 quadrature* approach, 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 sqrt(4 - x2) b/w limits of zero and two and you have a quarter of a circle of radius two ie. the area is going to be PI. This time I use the trapezoidal rule for generating sums, which conceptually looks like this :
Here 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 Simpson's Rule which uses a quadratic/parabolic approximation to the curve for each ordinate division, in detail implementing the so-called 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 vs iterations curve. See how much pain for how much gain ...... )
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
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 vs the exponent of the number of ordinate divisions :
... 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.14159265358979311600
it 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 vpaintegral command. Aha! A shiny new function to play with .....
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
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 + x2) 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 + x2) 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 + x2) = 1 - x2 + x4 - x6 + x8 + .... + (-1)k x2k + .....
.... 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 x2k integrates to x2k+1/(2k+1) } giving :
Integral[(1 + x2)-1] = x - x3/3 + x5/5 - x7/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't explain why the ratio of a circle's circumference to it's diameter is involved though ! What geometric construct, involving circles, yields (1 + x2)-1 as a solution ? Roger Penrose in his epic book The Road To Reality { chapter four Magical Complex Numbers } has an interesting take on all this, as he examines (1 + z2)-1 in the complex plane rather than merely (1 + x2)-1 on the real line.
Cheers, Mike.
( edit ) OK. Try this :
.... an actual tangent line to the circle whose radius is one, where (1 + x2)-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 :
Sn = 1 + 1/2 + 1/3 + 1/4 + ..... + 1/n
does 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 :
Sn = 1 + 1/4 + 1/9 + 1/16 + .... + 1/n2
does converge to the value of PI2/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 Grant Sanderson'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 :
"Answer = 1.64493405783457502523 ( 100000000 terms )"as it's best estimate on my machine before something breaks. Proving it in detail requires the use of complex numbers.
One could reasonably ask the question : for what value of p does
Sn = 1 + 1/2p + 1/3p + 1/4p + ..... + 1/np
such a series of partial sums converge ? Due to a method called the Comparison Test the series with 1/n3 as terms converges whereas 1/n0.5 doesn't. As above it doesn't for p= 1 and it does for p=2, so maybe some value of p in between is the limiting case ? It turns out that for any p > 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.
Next up there is the product by Mr John Wallis which converges to PI/2 as more multiplicative factors are included, and invented long before either Calculus or the Binomial Theorem were published :
PI/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 n tends to infinity :
((n!)2*(22n))/((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 ( e ) and various combinations like ePI - yes, that is e to the power of PI. But enough said for tonite methinks. :-)
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 thing. Look at
)
Now here's a thing. Look at these partial sums again :
Sn = 1 + 1/2p + 1/3p + 1/4p + ..... + 1/np
and now study in more detail the dependence of the final sum upon p, for p values really, really close to one but just above one. By the previous post these should always give convergent series. You get :
Where the data are plotted for successive runs of the summating machine. The first point on the bottom left is the PI2/6 limit. Then are plotted the ultimate sums after 100000000 terms each, for values of p closer and closer to 1 eg. the last point on the top right is for p = 1 + 10-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 19. 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 is 19 doing here ? It's a prime number ....
What I'm getting at is that one can define numbers by limiting processes.
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