Zeno’s Enduring Example
2500 years after settling a controversy and now as a reference frame origin, is the geometric series a key to improving math education?
Introduction
2500 years ago, Greek mathematicians had a problem with walking from one place to another. Physically, they were able to walk as well as we do today, perhaps better. Logically, however, they believed that every infinitely long list of numbers greater than zero sums to infinity. Therefore, it was a paradox when Zeno of Elea (circa 500 B.C.) pointed out that in order to walk from one place to another, you first have to walk half the distance, and then you have to walk half the remaining distance, and then you have to walk half of that remaining distance, and you continue halving the remaining distances an infinite number of times because no matter how small the remaining distance is … you still have to walk the first half of it. In his example, Zeno of Elea transformed even a short walking distance into an infinitely long list of halved remaining distances, all of which are greater than zero. And that was the problem: how can a distance be short when measured directly and also infinite when summed over its infinite list of halved remainders? The paradox revealed something was wrong with the assumption, widely held among Greek mathematicians at that time, that every infinitely long list of numbers greater than zero sums to infinity.
The series used in Zeno of Elea’s paradox, 1/2 + 1/4 + 1/8 + 1/16 + 1/32 + … , is now called a geometric series, as is every series of the form a + ar + ar² + ar³ + ar⁴ + … , where a is called the coefficient which is common to every term in the series and r is the common ratio that multiplies any term to generate the next term in the series. It is often a convenience in notation to set the series equal to the sum S and work with the geometric series
S = a + ar + ar² + ar³ + ar⁴ + … in its normalized form
S/a = 1 + r + r² + r³ + r⁴ + … or in its normalized vector form
S/a = [1 1 1 1 1 …][1 r r² r³ r⁴ …]ᵀ or in its normalized partial series form
Sₙ/a = 1 + r + r² + r³ + … + rⁿ, where n is the power (or degree) of the last term included in the partial sum Sₙ.
In this essay I propose that the same geometric series that resolved a controversy among Greek mathematicians 2500 years ago is now a key to improving today’s mathematics education. According to Jeff Hawkins on page 88 of his book, A Thousand Brains − a new theory of intelligence, “Being an expert is mostly about finding a good reference frame to arrange facts and observations.” In this essay, the geometric series is the origin of a proposed reference frame for arranging facts and observations about the body of knowledge of mathematics.
The phrase “body of knowledge” according to Wikipedia means “the complete set of concepts, terms and activities that make up a professional domain.” Given that the body of knowledge of mathematics covers so many concepts, it may be difficult to imagine “a good reference frame to arrange facts and observations” about the “complete set of concepts, terms and activities” of mathematics. In this essay I propose using a multi-floor building’s floor plan as an analogy. Specifically, I propose a floor plan for what is formally called the polynomial vector space (i.e., the set of all possible polynomials) and then I attempt to make the case that this polynomial vector space should be the ground floor entry point to all the other floors of mathematics.
Although this essay focuses on the polynomial vector space as the ground floor of the body of knowledge of mathematics, there are certainly many other floors, each with its own metric for measuring distance, in this building analogy. For example, perhaps a good second floor would hold all the concepts related to probability and statistics. And perhaps category theory could be used to reveal similarities between some of the many different floors. However, the focus of this essay is just the polynomial vector space ground floor as shown in the above figure where the origin of the log log plot is the geometric series and distances to other polynomial vector space landmarks (e.g., the complex Fourier series) are measured as the degrees of freedom away from the geometric series.
The hypothesis is that providing a good reference frame of the body of knowledge of mathematics would significantly increase everyone’s understanding of mathematics. Of course, a hypothesis needs to be tested to make it more than a supposition. Therefore, this essay’s conclusion proposes a way to test the hypothesis. Prior to that conclusion, this essay locates about a dozen landmarks in the proposed reference frame, occasionally referring to findings in the fields of neuroscience and psychology. For clarity, those references to neuroscience and psychology are in italics and highlighted by a vertical bar in the left margin, as follows.
Jeff Hawkins writes on page 50 of his book A Thousand Brains − a new theory of intelligence, “… we need to think of the neocortex as primarily processing reference frames. Most of the circuitry is there to create reference frames and track locations. Sensory input is of course essential. As I will explain in coming chapters, the brain builds models of the world by associating sensory input with locations in reference frames.
Why are reference frames so important? What does the brain gain from having them? First, a reference frame allows the brain to learn the structure of something. A coffee cup is a thing because it is composed of a set of features and surfaces arranged relative to each other in space. Similarly, a face is a nose, eyes, and a mouth arranged in relative positions. You need a reference frame to specify the relative positions and structure of objects.
Second, by defining an object using a reference frame, the brain can manipulate the entire object at once. For example, a car has many features arranged relative to each other. Once we learn a car, we can imaging what it looks like from different points of view or if were stretched in one dimension. To accomplish these feats, the brain only has to rotate or stretch the reference frame and all the features of the car rotate and stretch with it.
Third, a reference frame is needed to plan and create movements. Say my finger is touching the front of my phone and I want to press the power button at the top. If my brain knows the current location of my finger and the location of the power button, then it can calculate the movement needed to get my finger from its current location to the desired new one. A reference frame relative to the phone is needed to make this calculation.”
In addition to representing the reference frame origin in the ground floor of the body of knowledge of mathematics, the geometric series can serve as a consistent introduction to many different layers of abstraction in mathematics and those abstraction layers are this essay’s organization:
- quantities,
- numbers,
- algebra,
- geometry,
- functions,
- vector spaces, and
- matrix exponentials.
For brevity, this essay does not illustrate more than seven layers of abstraction but it could. Also for brevity, the software used in this essay is optional reading in the appendix. However, there is a lot of truth in the Latin proverb Docendo discimus meaning “by teaching, we learn” and software is a form of “teaching” a computer to do a task that frequently exposes the teacher’s many initially overlooked misunderstandings. Therefore, quickly browsing through the appendix would probably be a good exercise. Expanding upon that code in the appendix to “teach” your own computer something new would be even better.
To be clear, I am not a mathematician, a psychologist, a teacher, or a neuroscientist. However as a software engineer, I have often wondered why so many people (myself included) find mathematics more difficult to learn than other subjects. If the relative difficulty of learning math is due to its lack of a good reference frame and/or its many layers of abstraction, perhaps this essay’s proposed floor plan analogy and its approach of using the geometric series as a consistent and familiar foothold in the introduction to each abstraction layer may be a helpful structure for learning math in general, especially for high school students in the midst of forming their mental models about math.
Ideally, this essay will entice some actual mathematicians, psychologists, teachers, and neuroscientists into testing this essay’s hypothesis (i.e., that the consistent use of a geometric series as the reference frame origin when introducing math concepts to students would significantly improve their understanding of those topics) is true ... or not. However, there is historical evidence that consistent repetition works.
Wikipedia describes the illusory truth effect with these examples: “Although the truth effect has been demonstrated scientifically only in recent years, it is a phenomenon with which people have been familiar for millennia. One study notes that the Roman statesman Cato closed each of his speeches with a call to destroy Carthage, knowing that the repetition would breed agreement, and that Napoleon reportedly “said that there is only one figure in rhetoric of serious importance, namely, repetition”, whereby a repeated affirmation fixes itself in the mind “in such a way that it is accepted in the end as a demonstrated truth”. Others who have taken advantage of the truth effect have included Quintilian, Ronald Reagan, Bill Clinton, George W. Bush, Donald Trump, and Marcus Antonius in Shakespeare’s Julius Caesar. Advertising that repeats unfounded claims about a product may boost sales because some viewers may come to think that they heard the claims from an objective source. The truth effect is also used in news media and is a staple of political propaganda.”
Given that we are all biased to trust the familiar (unfortunately even when that familiarity is due to the repetition of false statements), why not employ that same bias to encourage confidence in math students? Consistently relating currently disparate math concepts back to the same geometric series would increase familiarity with unfamiliar math topics.
abstraction layer 1: quantities
In geometry, two triangles are similar when the three angles of one triangle are the same as the three angles in the other triangle. That happens by scaling the three side lengths of one triangle by some factor w to calculate the side lengths of the other triangle. Because the area of a triangle is A = bh/2, where A, b, and h are the triangle’s area, base, and height respectively, the area of a triangle with the side lengths scaled by w is A = (wb)(wh)/2 = w²bh/2.
In this essay, each term of a geometric series is represented as the area of a similar triangle. Therefore, in this essay, a sequence of similar triangles all have the same three angles and have areas that progressively change by the common ratio r. To account for that progressive change in similar triangle area by a factor of common ratio r, the similar triangle side lengths must progressively change by a factor of √r, because rA = rbh/2 = (wb)(wh)/2 = w²bh/2, or w=√r.
Although it is not unusual to represent quantities as areas (e.g., slices of a pie chart), it is unusual to represent quantities as partially overlapped (or occluded) areas. For example in Figure 1, the area of the large red triangle is 1 (the value of the first term in the normalized geometric series) because the triangle base is b=2, the triangle height is h=1, and the triangle area is A=bh/2 = 2/2 =1, but only the rightmost 1−r of that area is visible. Similarly, the area of the next largest (green) triangle is r (the value of the second term in the geometric series) because the triangle base is 2√r, the triangle height is √r, and the triangle area is A = bh/2 = 2√r√r/2=r, but only the rightmost r−r² of it is visible. Continuing, the area of the next largest (blue) triangle is r² (the value of the third term in the geometric series) because the triangle base is 2r, the triangle height is r, and the triangle area is A = bh/2 = 2rr/2=r², but only the rightmost r²−r³ of it is visible. From largest to smallest, the area of each successive geometric similar triangle decreases by a factor of r because the triangle base and height both decrease by a factor of √r. Also, the overlapping of area rⁿ by area rⁿ⁺¹ leaves only the rightmost rⁿ−rⁿ⁺¹ visible for all n ≥ 0.
This essay’s use of overlapped similar triangle areas to represent values of geometric series terms may at first seem awkward and unnecessary. However, that overlapping of similar triangles is the basis to revealing several geometric insights about the geometric series. In fact, as a sneak preview, the remainder of this essay is basically about manipulating these overlapped similar triangles:
- changing the shape of those triangle areas (i.e., removing the overlap),
- changing the orientation of those areas (i.e., rotating about the base),
- changing the height of those areas (i.e., varying the coefficients), and
- offering new perspectives about the sum of those areas.
According to the Scott Johnson’s Advanced Review in UCLA Baby Lab’s Development of visual perception, “… neonates are unable to perceive occlusion, and that occlusion perception emerges over the first several postnatal months. Interestingly, all of these effects in infants depend on the occluded stimulus moving behind the occluder, unlike adults who can perceive occlusion even with static images.” Dr. Johnson goes on to write in a Developmental Psychology article that “… perceptual completion in static occlusion displays has not been observed until 6–8 months.”
Although this essay’s representation of terms of the geometric series as partially occluded areas is unusual, the problem of recognizing partially occluded objects in the real world is important enough to warrant the early development (e.g., 6–8 months) of that capability.
abstraction layer 2: numbers
Zeno of Elea’s geometric series with coefficient a=1/2 and common ratio r=1/2 is the foundation of binary encoded approximations of fractions in digital computers. Concretely, the geometric series written in its normalized vector form is S/a = [1 1 1 1 1 …][1 r r² r³ r⁴ …]ᵀ. Keeping the column vector of basis functions [1 r r² r³ r⁴ …]ᵀ the same but generalizing the row vector
[1 1 1 1 1 …] so that each entry can be either a 0 or a 1 allows for an approximate encoding of any fraction. For example, the value v = 0.34375 is encoded as
v/a = [0 1 0 1 1 0 …][1 r r² r³ r⁴ r⁵ …]ᵀ where coefficient a = 1/2 and common ratio r = 1/2. Typically, the row vector is written in the more compact binary form v = 0.010110 which is 0.34375 in decimal.
Similarly, the partial geometric series with coefficient a=1 and common ratio r=2 is the foundation for binary encoded integers in digital computers. Again, the geometric series written in its normalized vector form is S/a =
[1 1 1 1 1 …][1 r r² r³ r⁴ …]ᵀ. Keeping the column vector of basis functions
[1 r r² r³ r⁴ …]ᵀ the same but generalizing the row vector [1 1 1 1 1 …] so that each entry can be either a 0 or a 1 allows for an encoding of any integer. For example, the value v = 151 is encoded as
v/a = [1 1 1 0 1 0 0 1 0 …][1 r r² r³ r⁴ r⁵ r⁶ r⁷ r⁸ …]ᵀ where coefficient a = 1 and common ratio r = 2. Typically, the row vector is written in reverse order (so that the most significant bit is first) in the more compact binary form v = …010010111 = 10010111 which is 151 in decimal.
As shown in the adjacent figure, the standard binary encoding of a 32-bit floating point number is a combination of a binary encoded integer and a binary encoded fraction, beginning at the most significant bit with
- the sign bit, followed by
- an 8-bit integer exponent field with an assumed offset of 127 (so a value of 127 represents an exponent value of 0) and with a base of 2 meaning that the exponent value specifies a bit shift of the fraction field, followed by
- a 23-bit fraction field with an assumed but not encoded 1 serving as the fraction’s most significant nonzero bit which would be in bit position 23 if it were encoded.
Building upon the previous example of 0.34375 having binary encoding of 0.010110, a floating point encoding (according to the IEEE 754 standard) of 0.34375 is
- the sign bit which is 0 because the number is not negative,
- an 8-bit integer exponent field which must specify a shift that counters the 2 bit left shift to get the original binary encoding from 0.010110 to 1.0110, and that counter shift to recover the original binary encoding is a right shift of 2 bits which is specified by an exponent value of 125 (because 125 − 127 = -2 which is a right shift of 2 bits) which in binary is 0111 1101,
- a 23-bit fraction field: .0110 0000 0000 0000 0000 000.
Although encoding floating point numbers by hand like this is possible, letting a computer do it is easier and less error prone. The following Julia code confirms the hand calculated floating point encoding of the number 0.34375:
julia> bitstring(Float32(0.34375))
"00111110101100000000000000000000"Decimal fractions that have repeating patterns that continue forever (e.g., 0.3333… or 0.09090909… or 0.123412341234…) can only be approximated when encoded as floating point numbers, but they can always be defined exactly as the ratio of two integers and those integers can be calculated using the geometric series. The trick is to use the closed form of the geometric series which is S = a/(1 − r) when |r|<1. That closed form of the geometric series will be derived in this essay’s discussion of abstraction layer 3 (algebra) and again in its discussion of abstraction layer 4 (geometry). However even before those derivations, you can confirm that S/a = 1/(1 − r) = 1 + r + r² + r³ + r⁴ + … when |r|<1 by long division of 1 by (1 − r).
As an example of using the geometric series to determine the numerator and denominator, the rational number 0.3333… can be written as the geometric series S = 0.3 + 0.3(1/10) + 0.3(1/10)² + 0.3(1/10)³ + … where coefficient a = 0.3 and common ratio r = 1/10. Therefore, |r|<1 and S = a/(1 − r) = (3/10)/(1 − (1/10)) = (3/10)/(9/10) = 3/9 = 1/3 = 0.333… . Similarly, 0.09090909… can be written as the geometric series S = 0.09 + 0.09(1/100) + 0.09(1/100)² + 0.09(1/100)³ + … where coefficient a = 0.09 and common ratio r = 1/100 resulting in the closed form S = a/(1 − r) = (9/100)/(1 − (1/100)) = (9/100)/(99/100) = 9/99 = 1/11 = 0.09090909… .
In addition to being able to calculate the exact numerator and denominator of rational numbers, and being at the core of the binary encoding used by most computers, the geometric series has another trick: its terms can rotate at harmonic frequencies. As shown in the adjacent animation, the geometric series can be extended into the set of complex numbers by setting the common ratio r to a complex number r = |r|eⁱᶿ, where |r| is the magnitude of that complex number and θ is its orientation angle in the complex plane. Making the common ratio r complex extends Figure 1 from 2-dimensions to 3-dimensions by adding a z-axis of imaginary numbers, directed out of the screen in Figure 1 and directed up in Figure 2b, so that the yz-plane is a complex plane. The quantities shown in Figure 1 rotate about the x-axis in Figure 2b. What makes this rotation particularly useful is that the similar triangle areas representing the terms of the complex geometric series rotate at different speeds, specifically at integer multiples of some fundamental frequency. The sum of all the rotating terms of every complex geometric series always traces out a perfect circle, provided that |r|<1.
People contemplate a wide variety of topics from esoteric topics such as encoded numbers to the widespread curiosity about the mechanism for thought itself. In the human brain, the neocortex is the outermost and largest part of the cerebral cortex which is the outermost and largest part of the cerebrum which is the outermost and largest part of the brain. The brain’s evolution added newer parts on top of older parts and they all affect each other. By volume, the neocortex is about 70% of the human brain. If straightened out, the neocortex is only about 2.5 millimeters (0.10 inch) thick but it has many folds and creases to pack more of itself into the skull.
abstraction layer 3: algebra
In addition to the previously mentioned normalized form and the partial series form of the geometric series, the closed forms are very practical in that they drastically reduce the number of operations required to calculate the series. The widely known derivation of the closed form of the geometric series subtracts off many of the inner terms like this:
Sₙ/a = 1 + r + r² + r³ + … + rⁿ⁻¹ + rⁿ,
rSₙ/a = r + r² + r³ + … + rⁿ + rⁿ⁺¹,
Sₙ/a−rSₙ/a = 1− rⁿ⁺¹,
(1−r)Sₙ/a = 1− rⁿ⁺¹,
Sₙ/a = (1−rⁿ⁺¹) / (1−r).
When r<0, the inner terms of the partial series are cancelled by addition instead of subtraction. Substituting p=−r>0 (where p is for positive) helps keep track of the signs.
Sₙ/a = 1 − p + p² − p³ + … − pⁿ⁻¹ + pⁿ (when n is even),
pSₙ/a = p − p² + p³ −… + pⁿ⁻¹ − pⁿ + pⁿ⁺¹ (when n is even),
Sₙ/a+pSₙ/a = 1 + pⁿ⁺¹ (when n is even),
Sₙ/a = (1 + pⁿ⁺¹)/(1+p) (when n is even),
Sₙ/a = (1 − rⁿ⁺¹)/(1−r) (when n is even).
Similarly,
Sₙ/a = 1 − p + p² − p³ + … + pⁿ⁻¹ − pⁿ (when n is odd),
pSₙ/a = p − p² + p³ − … − pⁿ⁻¹ + pⁿ − pⁿ⁺¹ (when n is odd),
Sₙ/a+pSₙ/a = 1− pⁿ⁺¹ (when n is odd),
Sₙ/a = (1 − pⁿ⁺¹)/(1+p) (when n is odd),
Sₙ/a = (1 − rⁿ⁺¹)/(1−r) (when n is odd).
Conveniently, Sₙ/a = (1 − rⁿ⁺¹)/(1 − r) is the same result when
- r>0 for even or odd n,
- r<0 for even n, and
- r<0 for odd n.
To extend the partial sum Sₙ/a to the full sum S/a, note that rⁿ⁺¹approaches 0 as n approaches infinity if and only if |r|<1. Then, the closed form of the geometric series is
S/a = 1 / (1 − r), but only within the range |r|<1.
That closed form (and its proof) still applies when the common ratio r is a complex number.
Although the neocortex is about 70% of your brain by volume, the other 30% is critical and often overrules your neocortex. For example, your neocortex might convince you to hold your breath. However, after about one minute for most people, the autonomous respiratory function of your brainstem will assert itself and put to rest your neocortex’s idea of you holding your breath.
STEM teachers typically focus on training the student’s neocortex and overlook the ancient but critically important 30% of the brain that isn’t the neocortex. In contrast, advertisers often focus entirely on that ancient 30% of the brain, manipulating our biases and emotions and our fascination with the unusual or unexpected. As a result, most advertisers are much more engaging than most STEM teachers.
abstraction layer 4: geometry
Observing a math problem from a slightly different perspective occasionally magically simplifies the solution. For example, the geometric series 4/9 + 4/9² + 4/9³ + 4/9⁴ + … , with coefficient a = 4/9 and common ratio r = 1/9, has a sum that is not obvious from looking at the expanded form of the series. However, using the algebraically derived geometric series closed form formula when |r|<1, the sum of 4/9 + 4/9² + 4/9³ + 4/9⁴ + … is easy to calculate as
S = a / (1 − r) = (4/9)/(1–1/9) = (4/9)/(8/9) = 1/2. The adjacent figure represents the value of each term of that geometric series as the area of a purple square, illustrating a perspective that makes the sum of 4/9 + 4/9² + 4/9³ + 4/9⁴ + … even easier to calculate. Concretely, by observing that the unit square is composed of an infinite number of L-shaped areas, each with an equal amount of purple and yellow area (i.e., 4 purple squares and 4 yellow squares in each L-shaped area), the total of the infinite number of purple areas must be half the area of the unit square, which is 1/2(1) = 1/2. Although Figure 4a’s geometric perspective helps quickly reveal the sum of that specific geometric series with coefficient a = 4/9 and common ratio r = 1/9, a geometric perspective of the geometric series closed form formula itself would help with visualizing every geometric series that converges. For that, we leave the squares behind and we return to the similar triangle area representations of the terms of the geometric series.
Removing all the overlaps among the similar triangles shown in Figure 1 reveals a simple geometric interpretation of the geometric partial series approaching the geometric series closed form: as the number of terms in the partial series approaches infinity, the area representing the series sum approaches the area of a triangle. Concretely, each term can be represented as an area of a trapezoid and those trapezoids aggregate into a larger trapezoid that converges to a triangle as more terms are added. There are several rewards for learning this geometric interpretation of the closed form of the geometric series:
- A trapezoid transforming into a triangle is a refreshing geometric interpretation of the well-known algebraic interpretation of the geometric series convergence.
- It is also nice to be able to see the shrinking triangle tip at the end of the aggregated trapezoid as being the partial geometric series approximation error (i.e., S−Sₙ, which is the quantity that the partial series Sₙ falls short of the full series S).
- Being able to geometrically interpret the well-known closed form of the geometric series is like finally being able to associate a face with a person’s name.
In the adjacent figure, the common ratio is in the range 0<r<1. The geometric derivation of the closed form of this geometric series has three main steps:
- (top plot) Represent the first n+1 terms of Sₙ/a as areas of overlapped similar right triangles as done in Figure 1.
- (middle plot) In the order from largest to smallest, remove each triangle’s overlapped area, which is always a fraction r of its area, and scale the remaining 1−r of the triangle’s non-overlapped area by 1/(1−r) so the area of the formerly overlapped triangle, now the area of a non-overlapped trapezoid, remains the same.
- (bottom plot) Aggregate the resulting n+1 non-overlapped trapezoids into a single non-overlapped trapezoid and calculate its area. The area of that aggregated trapezoid represents the value of the partial series. That area is equal to the outermost triangle minus the empty triangle tip:
Sₙ/a = (1−rⁿ⁺¹) / (1−r), which simplifies to S/a = 1/(1−r) when n approaches infinity and |r|<1.
When the common ratio is in the range −1<r<0, the similar triangle areas alternate between positive and negative areas as shown in the adjacent figure. Aggregating all the positive areas results in a geometric series with common ratio r². Similarly, aggregating all the negative areas also results in a geometric series with common ratio r² but the series is also scaled by r. Applying the S/a = 1/(1−r) derived for the case of 0<r<1 to these two geometric series with common ratio r² results in a total area of S/a = 1/(1−r²) + r/(1−r²) = (1+r)/(1−r²) = 1/(1−r) for the case of −1<r<0, which matches the result for the case of 0<r<1. Therefore the geometric derivation confirms the algebraic derivation of the closed form of the geometric series: S/a = 1/(1−r) when |r|<1. As an aside, in the geometric derivation when|r|≥1 the aggregate trapezoid representing the partial series no longer converges to a fixed area triangle.
As shown in the adjacent figure, the terms of a geometric series represented as areas of overlapped similar triangles or, equivalently, areas of non-overlapped trapezoids have an infinitely repeating self-similarity when magnified. In other words, continuously zooming into the areas representing the terms of the geometric series reveals that the smaller terms of the converging geometric series have the same shape (but different scale) as the larger terms.
With the previous geometric proof that the areas of overlapped similar triangles (e.g., Figure 1) can be transformed into the areas of non-overlapped trapezoids (e.g., Figure 4b), the rotating similar triangles of Figure 2b can be transformed into the rotating trapezoids in the adjacent Figure 4d. Note that the height of the trapezoids in Figure 4d are scaled by 0.5 so that the rotating trapezoids fit within the same plot limits of Figure 2b. Therefore, the heights of the rotating trapezoids in Figure 4e are actually twice the height of the rotating triangles in Figure 2b. Of course, the areas of the rotating triangles and the rotating trapezoids are the same because they both represent the same terms of the complex geometric series.
As more terms are added, the sum of the areas of the rotating triangles in Figure 2b as well as the sum of the areas of the rotating trapezoids in Figure 4e traces out a shape that gets closer to the shape of a perfect circle. Inversive geometry is an elegant approach to showing that every complex geometric series, S/a = 1/(1−r), traces a perfect circle when the number of terms approaches infinity, r is complex, and |r|<1.
In Figure 4e above, according to the geometric series closed form formula, S/a = 1/(1− eⁱᶿ/2). The denominator (i.e., 1 − eⁱᶿ/2) clearly traces a circle (shown as the green circle in the adjacent figure) of radius 1/2 that is centered at 1+i0 as θ varies from 0 to 2π. Inversive geometry proves that the inversion of that (green) circle is also a circle (the blue circle in the adjacent figure). In this case, the inversion circle is the unit circle (the red circle in the adjacent figure). Informally, the inversion flips the contents of the (red) inversion circle inside out,
- mapping points outside the red circle (e.g., point P) to points inside the red circle (e.g., point P’ short for point P inverse),
- mapping points inside the red circle (e.g., point Q) to points outside the red circle (e.g., point Q’ short for point Q inverse), and
- mapping points on the red circle (e.g., point R) to themselves,
where the mapping is constrained by |OP||OP’| = |OQ||OQ’| = |OR||OR’| = 1. Note that this constraint is for the normalized form of the inversion circle. Generally, the constraint is written as |OP||OP’| = |OQ||OQ’| = |OR||OR’| = rᵢₙᵥ² where rᵢₙᵥ is the radius of the inversion circle. Also note that the mapping constraint for inversion in a circle is a form of geometric mean, and can be expressed as part of a geometric sequence. For example, referring to Figure 4f, the increasing sequence |OQ|, rᵢₙᵥ, |OQ’| can be written as rᵢₙᵥ/α, rᵢₙᵥ, αrᵢₙᵥ where α>1 and equals the ratio |OQ’|/rᵢₙᵥ and is the common ratio in that partial sequence. For brevity this essay does not cover the theorems and proofs of inversive geometry. In general those theorems and proofs, which are thoroughly covered by this paper by Kozai and Libeskind, have recurring themes of similar triangles and inscribed angles.
This essay in general emphasizes a visual geometric interpretation of the geometric series even though there may be people who prefer an algebraic interpretation of the geometric series and there may be people with less visual acuity. Jeff Hawkins in his recent book A Thousand Brains − a new theory of intelligence on page 157, “The size of the neocortical regions varies considerably between people. For example,, region V1, the primary visual region, can be twice as big in some people as in others. V1 is the same thickness for everyone, but the area, and hence the number of columns, can vary. A person with a relatively small V1 and a person with a relatively large V1 both have normal vision and neither person is aware of the difference. There is a difference, however; a person with a large V1 has higher acuity, meaning they can see small things. This might be useful if you are a watchmaker, for example.”
Even for people with little or no sight, the concept of overlapped similar triangles transforming into non-overlapped trapezoids can be a helpful perspective when learning about geometric series.
Now, having covered encoded numbers and the geometric series and the complex geometric series, it is possible to start filling in this essay’s proposed reference frame of the body of knowledge of mathematics as shown in the adjacent figure. Note that a closed form is a constrained version of the series. Concretely, the closed forms of both the geometric series and the complex geometric series require that |r|<1 and the figure indicates this constraint by plotting the triangles representing the closed forms inside the circles representing the series. The figure’s legend describes the four items plotted. Along with their plotted locations, those four items are conceptually similar to real world landmarks that our brains model through the use of place cells in the brain’s hippocampus. Unlike real world distances that our brains model through the use of grid cells in the brain’s entorhinal cortex, the “distances” in Figure 4g measure degrees of freedom. For example, the common ratio of the geometric series is separated from the common ratio of the complex geometric series by a “distance” of one degree of freedom because the common ratio of the geometric series can be any number from the real number line (i.e., one degree of freedom or D.O.F.) and the common ratio of the complex geometric series can be a combination of any number from the real number line and any number from the imaginary number line (i.e., two degrees of freedom). How does the brain model a new type of distance?
In A Thousand Brains − a new theory of intelligence on page 62, Jeff Hawkins writes, “… to learn a complete model of something, you need both grid cells and place cells. Grid cells create a reference frame to specify locations and plan movements. But you also need sensed information represented by place cells, to associate sensory input with locations in the reference frame.
The mapping mechanisms in the neocortex are not an exact copy of ones in the old brain. Evidence suggests that the neocortex uses the same basic neural mechanisms, but it is different in several ways. It is as if nature stripped down the hippocampus and entorhinal cortex to a minimal form, made tens of thousands of copies, and arranged them side by side in cortical columns. That becoame the neocortex.”
abstraction layer 5: functions
As a quick review of the key points covered in this essay so far,
- the geometric series a + ar + ar² + ar³ + ar⁴ + … is fully specified by just two parameters: coefficient a and common ratio r,
- for a geometric perspective, the terms of the geometric series can be drawn as areas of overlapped similar triangles,
- those overlapped similar triangles can be transformed into non-overlapped trapezoids that, in aggregate, converge to a triangle as the number of geometric series terms increases,
- making the common ratio r a complex number rotates the terms of the complex geometric series at harmonic frequencies,
- those rotating terms can be drawn as overlapped similar triangles or non-overlapped trapezoids rotating about their bases, and
- provided |r|<1, the sum of those directed and rotated areas representing the terms of the complex geometric series traces a perfect circle as the angle (in the complex plane) of r traverses from 0 to 360 degrees.
Through all those extensions to the geometric series, coefficient a has not been extended ... until now. Allowing the coefficient of each term of the geometric series to vary independent of the other coefficients turns the geometric series into a power series and can be written as a₀ + a₁r + a₂r² + a₃r³ + a₄r⁴ + …
Geometrically, the terms of the power series can still be thought of as the overlapped similar triangles or the non-overlapped trapezoids of the geometric series except each term’s area is now scaled in height by its own coefficient. For the case of the complex geometric series, being able to vary the term coefficients allows the complex power series to trace more shapes than perfect circles. Note that if all the terms were rotating at the same speed, their sum would still trace a perfect circle even with different coefficients but, because the terms are rotating at integer multiples of some fundamental frequency, the complex power series can trace more than circles but not any arbitrary shapes. Informally, the shapes that can be traced by a complex power series are puffy and cloud-like, including circles of course for the case when all the coefficients are the same and the power series once again becomes the geometric series.
The question of what needs to be extended so that the complex power series can trace any arbitrary shape can be distilled down to the question of what needs to be extended so that the complex power series can trace a straight line segment. The answer to that is revealed in Euler’s formula:
eⁱᶿ = cosθ + i sinθ,
e⁻ⁱᶿ = cosθ − i sinθ,
eⁱᶿ + e⁻ⁱᶿ = 2cosθ,
cosθ = (eⁱᶿ + e⁻ⁱᶿ)/2.
From a geometry perspective, Euler’s formula shows that the sum of two identical vectors rotating in opposite directions traces a line segment. Basically, the counter rotation cancels a dimension of the movement in the complex plane. Removing the constraint that all terms must rotate in the same direction and removing the constraint that all terms must start rotating at the same orientation (or phase) enables the complex power series to extend to the complex Fourier series which can trace any arbitrary two dimensional closed figure (i.e., a figure that begins and ends at the same point).
Algebraically, the complex power series a₀ + a₁r + a₂r² + a₃r³ + a₄r⁴ + … is extended to the complex Fourier series f(θ, C) = … c₋₄e⁻ⁱ⁴ᶿ + c₋₃e⁻ⁱ³ᶿ + c₋₂e⁻ⁱ²ᶿ + c₋₁e⁻ⁱᶿ + c₀ + c₁eⁱᶿ + c₂eⁱ²ᶿ + c₃eⁱ³ᶿ + c₄eⁱ⁴ᶿ + … by allowing terms to rotate in either direction and by allowing the coefficients to be complex numbers so that they can specify the starting orientation, or phase, of each term’s rotation. The sum of those terms can be written as f(θ, C) to indicate a function f that maps to some single value depending on the rotation variable θ and the variable C which is the vector of possibly complex coefficients. Because f(θ, C) repeats itself in each 360 degree rotation of the variable θ, f(θ, C) is a periodic function.
That complex Fourier series is written in a slightly different notation in the top equation of the adjacent figure. First, the expanded form of the complex Fourier series f(θ, C) = … c₋₄e⁻ⁱ⁴ᶿ + c₋₃e⁻ⁱ³ᶿ + c₋₂e⁻ⁱ²ᶿ + c₋₁e⁻ⁱᶿ + c₀ + c₁eⁱᶿ + c₂eⁱ²ᶿ + c₃eⁱ³ᶿ + c₄eⁱ⁴ᶿ + … is shortened to the sigma summation notation. Second, the function’s input is t instead of θ, but θ is equal to 2πf₀t and that scaled input makes no difference to the function’s output. Finally, to slightly simplify the parameterization of the function and its integration, the period T is arbitrarily set to 1 which means f₀ (=1/T) is also 1. The bottom equation of Figure 5a can be confirmed to be true by substituting the expanded form of the complex Fourier series as f(t) into the integral and noting that all the rotating terms integrate to zero leaving just the integral of the constant cₙ from time 0 to 1, which is cₙ. A nice illustration of the rotating terms integrating to zero is in this 3blue1brown video.
As an example of the ability of the complex Fourier series to trace any two dimensional closed shape, the adjacent animation shows a complex Fourier series converging to a drawing of the letter ‘e’ (for exponential). All the terms rotating in the counter clockwise direction, which is the positive direction according to the right hand rule, are aggregated into a right “arm” (extending the right “hand” rule anatomy metaphor). Similarly, all the terms rotating in the clockwise direction are aggregated into a left arm. The sum of the endpoint of the right arm and the endpoint of the left arm is the approximation of the complex Fourier series shown as the symbol ‘+’ in the animation. Halfway from the origin to that ‘+’ symbol point is also halfway between the endpoint of the right arm and the endpoint of the left arm. That double midpoint is drawn as the symbol ‘o’ and, like a spider wrapping a recent victim, the animation’s left and right arms twirl around the ‘o’.
Having two arms drawing the complex Fourier series in Figure 5b seems like a natural extension to having one arm drawing the complex geometric series in Figure 4e. As shown in the back plane of Figure 5b, an alternative approach to aggregating the many rotating vectors into two arms is aggregating them into a single arm in which the adjacent segments alternate between rotating clockwise and rotating counter clockwise in a zigzag pattern. The two arm animation offers a perspective that reveals an inherent order to the coordinated motion. For example, the parameterization of the animation in Figure 5b traverses the large inner arc in the counter clockwise direction and then traverses the large outer arc in the clockwise direction. While traversing the large inner arc in the counter clockwise direction, the left arm (responsible for clockwise rotations) stays parked. Similarly, while traversing the large outer arc in the clockwise direction, the right arm (responsible for counter clockwise rotations) stays parked.
Although the two arm animation of the complex Fourier series may seem like a natural extension of the one arm animation of the complex geometric series, the complex geometric series has a common ratio r and the complex Fourier series does not (i.e., there is no mention of a common ratio in the complex Fourier series equations in Figure 5a and the common ratio is what defines the shape of the rotating trapezoid areas in Figure 5b). However, as shown in the adjacent figure, it is possible to calculate “phantom common ratios” for the complex Fourier series by performing a least squared error fit of the geometric series model arⁿ to the complex Fourier series coefficient data, excluding the coefficient of the constant term which specifies the arbitrary positioning of the letter ‘e’. The reward for this extension is that it enables a single 3D plot to show the shared context of the spatial data (the tracing of the letter ‘e’) and the frequency data (the height of the trapezoids rotating at different speeds) and the phase data (the starting orientations of those rotating trapezoids).
Looking at all the intricately coordinated motion in Figure 5b, it can be surprising to learn that all the coordination is explained and defined by the two simple equations in Figure 5a. Similarly, there is a lot of intricate coordination going on in your brain: your brain has approximately 150,000 cortical columns, each consisting of several hundred “minicolumns” each of which contains a little more than one hundred neurons each of which has thousands and sometimes tens of thousands of synapses to other neurons. However for the case of the human brain, we have not yet found simple principles that explain and define that intricate coordination.
According to Jeff Hawkins on page 39 of his book, A Thousand Brains, “People often say the brain is the most complicated thing in the universe. They conclude from this that there will not be a simple explanation for how it works, or that perhaps we will never understand it. The history of scientific discovery suggests they are wrong. Major discoveries are almost always preceded by bewildering complex observations. With the correct theoretical framework, the complexity does not disappear, but it no longer seems confusing or daunting.”
Adding the complex Fourier series “landmark” to this essay’s proposed reference frame of the body of knowledge of mathematics, the adjacent figure shows the complex Fourier series is a degree of freedom “distance” of 4n+1 away from the complex geometric series, where n is the highest power of a term in the partial complex Fourier series. The degree of freedom “distance” of 4n+1 is because the coefficient of the complex geometric series is a single real number chosen from the real number line (i.e., one degree of freedom in choice). In contrast, the complex Fourier series has 2n+1 terms: n terms rotating in the positive direction, n terms rotating in the negative direction, and one term not rotating at all. Being able to specify any initial phase of rotation for each term, each of those 2n+1 terms can have a complex coefficient resulting in 2*(2n+1) degrees of freedom in the choice of coefficients. The difference is 2*(2n+1) − 1 = 4n + 1.
abstraction layer 6: vector spaces
(TODO)
As shown in the adjacent figure of the normalized vector forms of the geometric series, the entries of the column vector are basis functions that match the terms of the normalized geometric sequence. Very much like a function, a row vector maps multiple input quantities (the column vector entries) into a single output quantity (the sum).
Removing the constraint that all the coefficients must be the same, the geometric series generalizes to the power series as shown in vector form in the adjacent figure, where the column vector variable is changed from r to x because the series is no longer a geometric series with a common ratio r. Any polynomial function can be specified with the correct linear combination of those basis functions in the column vector. In other words, the function space of the series shown in Figure 6b is the space of all polynomial functions. Some functions are in this function space even though they do not have the typical appearance of a polynomial function (e.g., 3x²+2x+7). For example, one of the ways to define the function eˣ is with the polynomial function eˣ = 1 + x/1! + x²/2! + x³/3! + x⁴/4! + … . However, many functions (e.g., sin(x)) are not in this function space and can only be approximated by a polynomial function. So there are potentially two dimensions of approximation: the first dimension is the partial series Sₙ approximating the full series S, and the second dimension is the polynomial function space approximating the actual function f(x). The Taylor series handles that second dimension of approximation by defining how to set the coefficients of a power series so that the polynomial function approximates the actual function f(x).
To calculate a function’s Taylor series coefficients, start by approximating the function with the expanded form of the polynomial f(x) ≈ c₀ + c₁x + c₂x² + c₃x³ + … , and set x=0 to determine c₀≈f(0). Then take the first derivative f⁽¹⁾(x) ≈ c₁ + 2c₂x + 3c₃x² + … , and set x=0 to determine c₁≈f⁽¹⁾(0). Then take the second derivative f⁽²⁾(x) ≈ 2c₂ + 3!c₃x + … , and set x=0 to determine c₂≈f⁽²⁾(0)/2. Then take the third derivative f⁽³⁾(x) ≈ 3!c₃ + 4!c₄x … , and set x=0 to determine c₃≈f⁽³⁾(0)/3!. Then take the fourth derivative f⁽⁴⁾(x) ≈ 4!c₄ + 5!c₅x … , and set x=0 to determine c₄≈f⁽⁴⁾(0)/4!. Provided the function f(x) is smooth enough for the derivatives to exist, continue in this manner to the k-th derivative to determine cₖ ≈ f⁽ᵏ⁾(0)/k! which is the Taylor series coefficient formula for approximating a section of f(x) that is centered at x=0.
Technically, those coefficients are called Maclaurin series coefficients because they are approximating a section of f(x) that is centered at x=0. Soon, that center point of approximation will be shifted to x₀ and the coefficients will then properly be called Taylor series coefficients. But before that, the following is an alternate approach to deriving those Maclaurin series coefficients. The alternate approach demonstrates that derivatives and integrals are linear operators (i.e., matrix multipliers) on polynomial functions (i.e., the coefficients in the row vector).
The top equation in the adjacent figure shows how a linear operator L multiplies the coefficients in the row vector. The next two equations are the derivative matrix D and the anti-derivative matrix D⁻¹. The linear operator L can be D or it can be D⁻¹ or it can be some combination of them or other linear operators. For example, L = D⁻¹D = I is the identity matrix. Of course, L = DD = D² is the second derivative and L = D⁻¹D⁻¹ = D⁻² is the second anti-derivative, more commonly called the second integral. Note that the generation of the indefinite integral constants is the reason that a one is prepended to the row vector and a zero is prepended to the column vector. Also note that at the risk of causing confusion, but following convention, the indefinite integral constants are represented by upper case ‘C’s and the power series coefficients are represented by lower case ‘c’s.
Similar to the previously shown derivation of the Maclaurin series coefficients that starts with a function being approximated by the expanded form of a power series, the following vector form derivation is also an exercise in calculating derivatives. Several derivative operators are shown in the adjacent image, the result of matrix multiplications. Referring to the row vector with the prepended 1 in the top equation of Figure 6c as R, and referring to that column vector with the prepended 0 as X, the starting point for the vector derivation of the Maclaurin series coefficients is the approximation f(x)≈RX evaluated at x=0 yields c₀≈f(0). The first derivative f⁽¹⁾(x)≈RDX evaluated at x=0 yields c₁≈f⁽¹⁾(0). The second derivative f⁽²⁾(x)≈RD²X evaluated at x=0 yields c₂≈f⁽²⁾(0)/2. The third derivative f⁽³⁾(x)≈RD³X evaluated at x=0 yields c₃≈f⁽³⁾(0)/6=f⁽³⁾(0)/3!. The fourth derivative f⁽⁴⁾(x)≈RD⁴X evaluated at x=0 yields c₄≈f⁽⁴⁾(0)/24=f⁽⁴⁾(0)/4!. Once again, provided the function f(x) is smooth enough for the derivatives to exist, continue in this manner to the k-th derivative f⁽ᵏ⁾(x)≈RDᵏX to determine cₖ≈f⁽ᵏ⁾(0)/k! when evaluated at x=0.
To approximate a section of f(x) that is instead centered at x=x₀, define a new function g(x) that is equivalent to f(x) but shifted along the x axis so that its origin is at f(x=x₀): g(x)=f(x+x₀). Repeating the now familiar exercise of calculating polynomial derivatives, but now for g(x), results in cₖ≈g⁽ᵏ⁾(0)/k! but g(0)=f(x₀) so cₖ≈f⁽ᵏ⁾(x₀)/k!. The polynomial variable in the Taylor series is, like the Maclaurin series, the distance from the center point of the approximation, but because of the shift of that center point to x=x₀, the distance from the center point for the Taylor series is x−x₀ instead of x. Therefore, the Taylor series approximation in vector form is the equation in Figure 6e.
(11 Dec 2021 note: The remainder of this essay is changing and I hope to finish the editing sometime this month … by the end of this year.)
Although the closed form of the geometric series (i.e., S/a = 1/(1−r)) can be derived with either algebra or geometry, the geometric derivation offers a better perspective on the rate of convergence varying with the common ratio r. The animation in the adjacent figure shows five normalized geometric series (common ratio r=−0.74, −0.37, 0, 0.37, 0.74) converging to the closed form as the number of terms in the partial series increases. As expected, the geometric series with common ratio r=0 requires only one term to converge. Also as expected, the trapezoids get thinner as |r| increases which increases the number of terms needed to converge. However given the asymmetry of S/a = 1/(1−r), the cyan curve in the animation, the perfect symmetry about r=0 of the rate of convergence can be surprising. Algebraically, (S−Sₙ)/S = rⁿ⁺¹ which is the same regardless of the sign of r when n is odd. One way to think of this is that near r=1 the terms (all with magnitudes less than 1) contribute to an infinite sum, and near r=−1 the terms (again all with magnitudes less than 1) contribute to an infinite slope. In the animation for larger values of n, Sₙ/a is too large to be in view near r=1, but near r=−1 Sₙ/a continues to be in view and you can see the little tail flipping back and forth getting closer and closer to an infinite slope.
Although a summation in the form of a row vector times a column vector (e.g., Figure 6a) implies a vector space, and a vector space with an infinite number of basis functions (e.g., the space of all polynomial functions) is a little more mysterious than a vector space with basis vectors of just a few numbers, the inherent mystery of vector spaces is due to the definition not specifying any particular basis at all. Instead, a vector space is required to satisfy eight axioms:
- vector addition is associative (e.g., u + (v + w) = (u + v) + w),
- vector addition is commutative (e.g., u + v = v + u),
- vector addition has identity element (e.g., v + 0 = v),
- vector addition has inverse elements (e.g., v + (−v) = 0),
- scalar and field multiplication compatible (e.g., a(bv) = (ab)v),
- scalar multiplication has identity element (e.g., 1v = v),
- scalar multiplication is distributive with respect to vector addition (e.g., a(u + v) = au + av),
- scalar multiplication is distributive with respect to field addition (e.g., (a + b)v = av + bv).
As an example of the vast range of possible vector spaces, I am currently attempting to design a vector space of categories as part of a new type of dictionary. It is still too soon to tell whether this attempt will be successful or not.
According to Jeff Hawkins on page 88 of his book, A Thousand Brains, “Being an expert is mostly about finding a good reference frame to arrange facts and observations. Albert Einstein started with the same facts as his contemporaries. However, he found a better way to arrange them, a better reference frame, that permitted him to see analogies and make predictions that were surprising. What is most fascinating about Einstein’s discoveries related to special relativity is that the reference frames he used to make them were everyday objects. He though about trains, people, and flashlights. He started with the empirical observations of scientists, such as the absolute speed of light, and used everyday reference frames to deduce the equations of special relativity. Because of this, almost anyone can follow his logic and understand how he made his discoveries. In contrast, Einstein’s general theory of relativity required reference frames based on mathematical concepts called field equations, which are not easily related to everyday objects. Einstein found this much harder to understand, as does pretty much everyone else.”
abstraction layer 7: matrix exponentials
(TODO)
conclusion
In a 1974 Caltech commencement address, Richard Feynman (a winner of the 1965 Nobel Prize in Physics) talked about cargo cult science, a phrase he used to describe a story about pseudoscience. “In the South Seas there is a Cargo Cult of people. During the war they saw airplanes land with lots of good materials, and they want the same thing to happen now. So they’ve arranged to make things like runways, to put fires along the sides of the runways, to make a wooden hut for a man to sit in, with two wooden pieces on his head like headphones and bars of bamboo sticking out like antennas — he’s the controller — and they wait for the airplanes to land. They’re doing everything right. The form is perfect. It looks exactly the way it looked before. But it doesn’t work. No airplanes land. So I call these things Cargo Cult Science, because they follow all the apparent precepts and forms of scientific investigation, but they’re missing something essential, because the planes don’t land.
[…]there is one feature I notice that is generally missing in Cargo Cult Science. […] It’s a kind of scientific integrity, a principle of scientific thought that corresponds to a kind of utter honesty — a kind of leaning over backwards. For example, if you’re doing an experiment, you should report everything that you think might make it invalid — not only what you think is right about it: other causes that could possibly explain your results; and things you thought of that you’ve eliminated by some other experiment, and how they worked — to make sure the other fellow can tell they have been eliminated.
Details that could throw doubt on your interpretation must be given, if you know them. You must do the best you can — if you know anything at all wrong, or possibly wrong — to explain it. If you make a theory, for example, and advertise it, or put it out, then you must also put down all the facts that disagree with it, as well as those that agree with it.”
“But then I began to think, what else is there that we believe? (And I thought then about the witch doctors, and how easy it would have been to check on them by noticing that nothing really worked.) So I found things that even more people believe, such as that we have some knowledge of how to educate. There are big schools of reading methods and mathematics methods, and so forth, but if you notice, you’ll see the reading scores keep going down — or hardly going up — in spite of the fact that we continually use these same people to improve the methods. […]
Yet these things are said to be scientific. We study them. And I think ordinary people with commonsense ideas are intimidated by this pseudoscience. A teacher who has some good idea of how to teach her children to read is forced by the school system to do it some other way — or is even fooled by the school system into thinking that her method is not necessarily a good one.”
And in this essay, I propose a method for making mathematics easier to learn. The method may or may not be a good one. Given that we don’t even yet understand the brain, every proposed education improvement should be met with skepticism. So how can this proposed method be tested to show that it actually works and is not another example of cargo cult science? I think there are two requirements to investigate it:
- Don’t interfere. Let the mathematics teachers in the investigation choose to use the proposed method or not. Supply a large poster of the proposed reference map of the body of knowledge of mathematics to be posted off to the side in the classroom, away from where the teacher usually addresses the class. Whenever the teacher stands next to the out-of-the-way poster, a smartphone will automatically record video that is transcribed.
- Account for the variability of people. Students have a wide range of abilities. So do teachers. To establish a statistical baseline of the abilities of the teachers and the students, only teachers with at least five years of experience teaching the same class at the same school should be included in the investigation. Then correlate changes in student performance on mathematics tests with the amount of references (in the transcribed video) to the proposed reference map. In other words, the proposed investigation would generate a dose response plot, like in drug investigations but here the dose is references to the proposed reference map and the response is scores on mathematics tests.
Appendices
Writing or exploring an algorithm that implements a math concept often leads to a better understanding of that concept. The following appendices list the source code that generates the figures and the animation frames in this essay. I initially used GNU Octave to generate the animation frames and that source code is in appendix A. However, to take advantage of the efficiency of its state-of-the-art LLVM compiler which turned out to generate the animation frames four times faster (e.g., in 30 minutes instead of 120 minutes), I ported the GNU Octave application to Julia and the Julia source code (tested using Julia version 1.7.0) for generating the frames of the complex Fourier series animation is in appendix B. The Julia source code for generating the frames of the complex geometric series animation is in appendix C. The little Julia code snippet that distributes the complex Fourier series animation frame generation to six CPU cores is in appendix D. (This parallel processing reduced the time to generate the animation frames to 11.1 minutes, a tenth of the time GNU Octave took to generate similar animation frames.) The Julia source code for generating miscellaneous plots (e.g., the geometric series self-similarity animation and the figure of the reference frame of the body of knowledge of mathematics) is in appendix E. The Julia source code for generating the animation of quadratic and cubic Bezier curves is in appendix F.
Appendix A. GNU Octave source code
The following GNU Octave code generates the frames for the animation of the complex Fourier series drawing of the letter ‘e’.
% define some plotting parameters
nT = 801 % number of Time samples
nTM = nT - 1; % number of Time sample Midpoints
T = linspace(0, 1, nT); % Time samples
nT4 = nTM/4;
dT = 1/nTM; % difference between Time samples
TM = T(1:end-1) + dT/2; % Time sample Midpoints
THETAM = 2*pi*TM; % angle samples
XZ = zeros(1,nT); % X Zeros
CT = ['r'; 'g'; 'b'; 'y'; 'w'; 'k']; % Color TabletrapScale = 0.5; % trapezoid scale
dCross = 0.15; % diameter of '+' cross-hair
yOffset = 0.5; % y offset of F(t)
nAng = 15;
ANGDOT = linspace(0,2*pi,nAng);
XDOT = zeros(1,nAng);
YDOT = dCross/pi*cos(ANGDOT);
ZDOT = dCross/pi*sin(ANGDOT);% define parameterized trapezoid
wChar = 0.1;
rInner = 0.5;
rOuter = rInner+wChar;
cy = 0.4;
thetaChar = 7*pi/4;
F = zeros(1, 2*nT); % left half real, right half imaginary
M = zeros(2,1); % slope (Re and Im)P = zeros(1,7); % Perimeter segments
P(1) = cy + rInner;
P(2) = rInner * thetaChar;
P(3) = wChar;
P(4) = rOuter * thetaChar;
P(5) = wChar;
P(6) = cy + rOuter;
P(7) = wChar;
pSum = sum(P);
TS = P/pSum; % Time Segments
TSC = cumsum(TS); % Time Segments CumulativeV = zeros(2, 7); % Vertices
V(1,1) = cy + rInner;
V(2,1) = 0;
V(1,2) = -rInner + rInner * cos(thetaChar);
V(2,2) = rInner * sin(thetaChar);
V(1,3) = wChar * cos(thetaChar);
V(2,3) = wChar * sin(thetaChar);
V(1,4) = -rOuter * cos(thetaChar) + rOuter;
V(2,4) = -rOuter * sin(thetaChar);
V(1,5) = 0;
V(2,5) = -wChar;
V(1,6) = -cy - rOuter;
V(2,6) = 0;
V(1,7) = 0;
V(2,7) = wChar;
VC = cumsum(V, 2); % Vertices CumulativeiSegment = 1;
i0 = 1;
tOffset = 0;
nTS = 1 + floor(TS(1)/dT);
i1 = i0 + nTS - 1;
M = V(:,1) ./ TS(1);
F(i0:i1) = M(1) * T(i0:i1);
F(nT+i0:nT+i1) = 0;++iSegment; % 2
i0 = i1 + 1;
tOffset = T(i0) - TSC(1);
nTS = 1 + floor((TS(2) - tOffset)/dT);
i1 = i0 + nTS - 1;
period = 8*TS(2)/7;
THETA = 2*pi/period*(T(i0:i1) - TSC(1));
F(i0:i1) = cy + rInner*cos(THETA);
F(nT+i0:nT+i1) = rInner*sin(THETA);++iSegment; % 3
i0 = i1 + 1;
tOffset = T(i0) - TSC(2);
nTS = 1 + floor((TS(3) - tOffset)/dT);
i1 = i0 + nTS - 1;
M = V(:,3) ./ TS(3);
F(i0:i1) = VC(1, 2) + M(1) * ...
(T(i0:i1) - TSC(2));
F(nT+i0:nT+i1) = VC(2, 2) + M(2) * ...
(T(i0:i1) - TSC(2));++iSegment; % 4
i0 = i1 + 1;
tOffset = T(i0) - TSC(3);
nTS = 1 + floor((TS(4) - tOffset)/dT);
i1 = i0 + nTS - 1;
period = 8*TS(4)/7;
THETA = thetaChar - 2*pi/period*(T(i0:i1) - TSC(3));
F(i0:i1) = cy + rOuter*cos(THETA);
F(nT+i0:nT+i1) = rOuter*sin(THETA);++iSegment; % 5
i0 = i1 + 1;
tOffset = T(i0) - TSC(4);
nTS = 1 + floor((TS(5) - tOffset)/dT);
i1 = i0 + nTS - 1;
M = V(:,5) ./ TS(5);
F(i0:i1) = VC(1, 4);
F(nT+i0:nT+i1) = VC(2, 4) + M(2) * ...
(T(i0:i1) - TSC(4));++iSegment; % 6
i0 = i1 + 1;
tOffset = T(i0) - TSC(5);
nTS = 1 + floor((TS(6) - tOffset)/dT);
i1 = i0 + nTS - 1;
M = V(:,6) ./ TS(6);
F(i0:i1) = VC(1, 5) + M(1) * ...
(T(i0:i1) - TSC(5));
F(nT+i0:nT+i1) = VC(2, 5);++iSegment; % 7
i0 = i1 + 1;
tOffset = T(i0) - TSC(6);
nTS = 1 + floor((TS(7) - tOffset)/dT);
i1 = i0 + nTS - 1;
M = V(:,7) ./ TS(7);
F(i0:i1) = VC(1, 6);
F(nT+i0:nT+i1) = VC(2, 6) + M(2) * ...
(T(i0:i1) - TSC(6));% apply y offset shift then take midpoints
F(1:nT) += yOffset;
FM = [ ...
(F(1:nT-1) + F(2:nT))/2 ...
(F(nT+1:end-1) + F(nT+2:end))/2 ];% integrate complex coefficient
NTERM = [1:8, 12:4:24];
nTerm = size(NTERM,2)
nTermMax = NTERM(end)
FEM = zeros(1, 2*nTM);
CR = zeros(5, 2*nTermMax+1); % coefficient Real component
CI = zeros(5, 2*nTermMax+1); % coefficient Imag component
iA = 0;
for iTerm=-nTermMax:nTermMax
iA++;
EM = [cos(-iTerm*THETAM) sin(-iTerm*THETAM)];
FEM(1:nTM) = ...
FM(1:nTM) .* EM(1:nTM) - ...
FM(nTM+1:end) .* EM(nTM+1:end);
FEM(nTM+1:end) = ...
FM(1:nTM) .* EM(nTM+1:end) + ...
FM(nTM+1:end) .* EM(1:nTM);
V = reshape(FEM .* dT, nTM, 2);
VR = reshape(V(:,1),nT4,4);
CR(1:4,iA) = (sum(VR))';
VI = reshape(V(:,2),nT4,4);
CI(1:4,iA) = (sum(VI))';
endfor;
CR(5,:) = sum(CR(1:4,:));
CI(5,:) = sum(CI(1:4,:));% calculate common ratio r that best fits coefficients
A = zeros(2, 1); % [a ccw terms; a cw terms]
R = zeros(2, 1); % [r ccw terms; r cw terms]
I_CMAG4 = zeros(2*nTermMax+1,5); % columns 3 thru 5 are for testing
I_CMAG4(:,1) = -nTermMax:nTermMax;
I_CMAG4(:,2) = sqrt(CR(5,:).^2 + CI(5,:).^2);
X = [ones(nTermMax,1) I_CMAG4(nTermMax+2:end,1)];
Y = log(I_CMAG4(nTermMax+2:end,2));
PHI = (X'*X) \ (X'*Y); % is faster than PHI = inv(X'*X)*X'*Y
EPHI = exp(PHI);
A(1) = EPHI(1);
R(1) = EPHI(2);
X(:,2) = abs(I_CMAG4(1:nTermMax,1));
Y = log(I_CMAG4(1:nTermMax,2));
PHI = (X'*X) \ (X'*Y); % is faster than PHI = inv(X'*X)*X'*Y
EPHI = exp(PHI);
A(2) = EPHI(1);
R(2) = EPHI(2);
Rclf;
figure(1);
plot(I_CMAG4(nTermMax+2:end,1), ...
log(A(1)*R(1).^I_CMAG4(nTermMax+2:end,1)), ...
'k-');
box off;
grid on;
hold on;
plot(I_CMAG4(1:nTermMax,1), ...
log(A(2)*R(2).^abs(I_CMAG4(1:nTermMax,1))), ...
'k-');
ICOLOR = mod(abs(I_CMAG4(:,1)),3) == 0;
plot(I_CMAG4(ICOLOR,1), log(I_CMAG4(ICOLOR,2)), ...
'ko', ...
'markerfacecolor', 'r', ...
'markersize', 5);
ICOLOR = mod(abs(I_CMAG4(:,1)),3) == 1;
plot(I_CMAG4(ICOLOR,1), log(I_CMAG4(ICOLOR,2)), ...
'ko', ...
'markerfacecolor', 'g', ...
'markersize', 5);
ICOLOR = mod(abs(I_CMAG4(:,1)),3) == 2;
plot(I_CMAG4(ICOLOR,1), log(I_CMAG4(ICOLOR,2)), ...
'ko', ...
'markerfacecolor', 'b', ...
'markersize', 5);
xlabel("coefficient index n");
ylabel("ln(coefficient magnitude)");
strTitle1 = "LSE fit to coefficient magnitude data";
strTitle2 = sprintf("+ index common ratio r = %.4f", R(1));
strTitle3 = sprintf("- index common ratio r = %.4f", R(2));
title({strTitle1, strTitle2, strTitle3});
hold off;
print("-dpng", "-r400", "-Farial:10", '_pco.png');% plot coefficient magnitudes as trapezoids
SMB = zeros(2, 1); % Sum M*B; dims 1: +/- coefficients
SB = zeros(2, 1); % Sum B; dims 1: +/- coefficients
TH = zeros(2, nTermMax+1); % Trapezoid Height; dims 1: +/- coefficients; dims2: trapezoid height (1-based coefficient index)
X = zeros(2, 2); % dims 1: +/- coefficients; dims 2: x of tall/short height
M = zeros(2, 1); % trapezoid slope; dims 1: +/- coefficients
B = zeros(2, 1); % trapezoid base width; dims 1: +/- coefficientsX(:,2) = 2 * (1 .- R .^ 0.5);
TH(:,1) = [I_CMAG4(nTermMax+1,2); I_CMAG4(nTermMax+1,2)] / 2; % split DC component into +/- coefficients
TH(:,1) = TH(:,1) ./ (1 .- R); % trapezoid height: dims 1: +/- coefficients
XTRAP = [ ...
X(1,1); X(1,1); X(1,2); X(1,2); X(1,1); ...
NaN; ...
-X(2,1); -X(2,1); -X(2,2);, -X(2,2); -X(2,1)];
YTRAP = [ ...
0; TH(1,1); TH(1,1)*R(1)^0.5; 0; 0; ...
NaN; ...
0; TH(2,1); TH(2,1)*R(2)^0.5; 0; 0];
I_CMAG4(nTermMax+1,3) = sum(TH(:,1) .* (1 .- R), 1);xlim([-2,2]);
set(gca, 'xtick', -2:2:2);
plot(XTRAP, YTRAP, 'r-');box off;
grid on;
hold on;
for jTerm=1:nTermMax
iColor = mod(jTerm,3) + 1; X(:,1) = X(:,2);
X(:,2) = 2 * (1 .- R .^ ((jTerm+1)/2));
B = abs(diff(X, 1, 2));
TH(:,jTerm+1) = [ ...
I_CMAG4(nTermMax+1+jTerm,2); ...
I_CMAG4(nTermMax+1-jTerm,2)] ./ ...
(R .^ (jTerm/2) .* (1 .- R));
M = TH(:,jTerm+1) .* (R .- 1) ./ B;
XTRAP = [ ...
X(1,1); X(1,1); X(1,2); X(1,2);, X(1,1);
NaN;
-X(2,1); -X(2,1); -X(2,2); -X(2,2); -X(2,1)];
YTRAP = [ ...
0; TH(1,jTerm+1); TH(1,jTerm+1)*R(1)^0.5; 0; 0;
NaN;
0; TH(2,jTerm+1); TH(2,jTerm+1)*R(2)^0.5; 0; 0];
plot(XTRAP, YTRAP, CT(iColor));
SB += B;
SMB += M .* B;
I_CMAG4(nTermMax+1+jTerm,3) = TH(1,jTerm+1) * R(1)^(jTerm/2) * (1 - R(1));
I_CMAG4(nTermMax+1-jTerm,3) = TH(2,jTerm+1) * R(2)^(jTerm/2) * (1 - R(2));
endMWA = SMB ./ SB % M (i.e., slope) Weighted Average
X = [0; 2; NaN; 0; -2];
Y = [-MWA(1)*2; 0; NaN; -MWA(2)*2; 0];
plot(X, Y, ...
'color', [1.0, 0.65, 0]);
hold off;
print("-dpng", "-r400", "-Farial:10", '_pto.png');
XMWA = [0, 2, NaN, 0, -2];
YMWA = [-trapScale*MWA(1)*2, 0, NaN, -trapScale*MWA(2)*2, 0];
ZMWA = [0, 0, NaN, 0, 0];% calculate partial series terms
PS = zeros(2*nTermMax+1, 2*nTM); % Partial Sum
iTerm0 = nTermMax + 1;
PS(1,1:nTM) = CR(5,iTerm0); % DC component doesn't rotate
PS(1,nTM+1:end) = CI(5,iTerm0); % DC component doesn't rotate
for jTerm=1:nTermMax
iPS = 2 * jTerm;
iTerm = iTerm0 + jTerm;
EM = [cos(jTerm*THETAM)' sin(jTerm*THETAM)']; % nTM x 2
PS(iPS,1:nTM) = ...
CR(5,iTerm) .* EM(:,1) - ...
CI(5,iTerm) .* EM(:,2);
PS(iPS,nTM+1:end) = ...
CI(5,iTerm) .* EM(:,1) + ...
CR(5,iTerm) .* EM(:,2); iPS = 2 * jTerm + 1;
iTerm = iTerm0 - jTerm;
EM = [cos(jTerm*THETAM)' sin(-jTerm*THETAM)']; % nTM x 2
PS(iPS,1:nTM) = ...
CR(5,iTerm) .* EM(:,1) - ...
CI(5,iTerm) .* EM(:,2);
PS(iPS,nTM+1:end) = ...
CI(5,iTerm) .* EM(:,1) + ...
CR(5,iTerm) .* EM(:,2);
endfor;
PSC = cumsum(PS, 1);% split terms according to direction of rotation
PSP = zeros(nTermMax+1, 2*nTM); % Partial Sum Positive rotation
PSP(1,:) = PS(1,:) / 2;
PSP(2:end,:) = PS(2:2:2*nTermMax+1,:);
PSPC = cumsum(PSP); % Partial Sum Positive rotation Cumulative
PSN = zeros(nTermMax+1, 2*nTM); % Partial Sum Negative rotation
PSN(1,:) = PS(1,:) / 2;
PSN(2:end,:) = PS(3:2:2*nTermMax+1,:);
PSNC = cumsum(PSN); % Partial Sum Negative rotation Cumulative
PST = PSP + PSN; % Partial Sum Total
PSTC = cumsum(PST, 1); % Partial Sum Total Cumulative% test PSP and PSN in column 4 of I_CMAG4
kTM = 50; % test should work for any value of kTM where 1<=kTM<=nTM
I_CMAG4(nTermMax+1,4) = (PSP(1,kTM)^2 + PSP(1,kTM+nTM)^2)^0.5;
I_CMAG4(nTermMax+1,4) += (PSN(1,kTM)^2 + PSN(1,kTM+nTM)^2)^0.5;
for jTerm=2:nTermMax+1
I_CMAG4(nTermMax+jTerm,4) = (PSP(jTerm,kTM)^2 + PSP(jTerm,kTM+nTM)^2)^0.5;
I_CMAG4(nTermMax+2-jTerm,4) = (PSN(jTerm,kTM)^2 + PSN(jTerm,kTM+nTM)^2)^0.5;
end% animate
figure(3, 'position', [600 400 800 450]);
FRAME0 = nTM/2 * ones(nTerm, 1);
FRAME0(1) = 0;
FRAME0(9) *= 2;
FRAME0 = cumsum(FRAME0, 1);
iTerm = 0;
for kTerm=NTERM
iTerm++;
iFrame = FRAME0(iTerm);
for iTM=1:nTM % occasionally skip every other frame for speed
if ((kTerm!=8) && (kTerm!=24))
if (mod(iTM,2)==0)
continue;
endif;
endif;
% setup plot figure
clf;
figure(3, 'position', [600 400 800 450]);
daspect([0.5 0.4 4]);
xlim([-2 2]);
ylim([-1 2]);
zlim([-1 1]);
set(gca, "fontsize", 12);
set(gca, 'xtick', -2:2:2);
set(gca, 'ytick', -1:2);
set(gca, 'ztick', -1:1);
grid on;
hold on;
xlabel('trapezoid baseline');
ylabel('Re {\bfS_n}');
zlabel('Im {\bfS_n}');
axis('equal');
view([72.5 20]);
strTitle1 = sprintf( ...
"complex Fourier series drawing of 'e'");
strTitle2 = sprintf( ...
"(trapezoid scale=%.2f, t=%.4f, n=+/-%d)", ...
trapScale, iTM/nTM, kTerm);
title({strTitle1,strTitle2}); % plot 'e' (closed form) in x=-2 plane
plot3(XZ-2, F(1:nT), F(nT+1:end), 'c-'); % plot weighted average trapezoid slopes (MWA)
plot3(XMWA, YMWA, ZMWA, ':', ...
'color', [1.0, 0.65, 0]); % plot zigzag aggregated vectors (x=-2 plane)
XVEC = [-2 -2];
YVEC = [0 PS(1,iTM)];
ZVEC = [0 PS(1,iTM+nTM)];
plot3(XVEC, YVEC, ZVEC, CT(1), ...
'linewidth', 1);
XVEC = [-2, -2, -2];
for iPS=2:2:2*kTerm
iColor = mod(floor(iPS/2),3) + 1;
YVEC = [PSC(iPS-1,iTM), PSC(iPS,iTM), PSC(iPS+1,iTM)];
ZVEC = [PSC(iPS-1,iTM+nTM), PSC(iPS,iTM+nTM), PSC(iPS+1,iTM+nTM)];
plot3(XVEC, YVEC, ZVEC, CT(iColor), ...
'linewidth', 1);
endfor; % plot axes
plot3(xlim(), [0 0], [0 0], '-', ...
'color', 0.85+[0 0 0]);
plot3([2 2], ylim(), [0 0], '-', ...
'color', 0.85+[0 0 0]);
plot3([2 2], [0 0], zlim(), '-', ...
'color', 0.85+[0 0 0]); % plot trapezoids
for iPS=1:kTerm+1
iColor = mod(iPS-1,3) + 1;
S = R .^ ((iPS-1)/2) .* (1 .- R); % plot trapezoids (mid planes)
XTRAP = [ ...
2*(1-R(1)^((iPS-1)/2));
2*(1-R(1)^((iPS-1)/2));
2*(1-R(1)^(iPS/2));
2*(1-R(1)^(iPS/2));
2*(1-R(1)^((iPS-1)/2));
NaN;
-2*(1-R(2)^((iPS-1)/2));
-2*(1-R(2)^((iPS-1)/2));
-2*(1-R(2)^(iPS/2));
-2*(1-R(2)^(iPS/2));
-2*(1-R(2)^((iPS-1)/2))];
YTRAP = [ ...
0;
trapScale*PSP(iPS,iTM)/S(1);
trapScale*PSP(iPS,iTM)/S(1) * R(1);
0;
0;
NaN;
0;
trapScale*PSN(iPS,iTM)/S(2);
trapScale*PSN(iPS,iTM)/S(2) * R(2);
0;
0];
ZTRAP = [ ...
0;
trapScale*PSP(iPS,iTM+nTM)/S(1);
trapScale*PSP(iPS,iTM+nTM)/S(1) * R(1);
0;
0;
NaN;
0;
trapScale*PSN(iPS,iTM+nTM)/S(2);
trapScale*PSN(iPS,iTM+nTM)/S(2) * R(2);
0;
0];
plot3(XTRAP, YTRAP, ZTRAP, CT(iColor), ...
'linewidth', 0.3);
% fill test column 5 of I_CMAG4
if ((kTerm == nTermMax) && (iTM == 1))
B = 2 .* R.^(iPS-1) .* (1 .- R.^0.5);
I = [2; 8];
H = ((YTRAP(I).^2 + ZTRAP(I).^2).^0.5) ./ trapScale;
if (iPS > 1)
I = [nTermMax+iPS; nTermMax+2-iPS];
I_CMAG4(I,5) = H .* R .^ ((iPS-1)/2) .* (1 .- R);
else
SDC = H .* B .* (1 .+ R.^0.5) ./ 2;
I_CMAG4(nTermMax+1,5) = sum(SDC);
end
end
endfor; % plot 'e' (closed form) in x=2 plane
plot3(XZ+2, F(1:nT), F(nT+1:end), 'c-'); % plot dotted trail of partial sum
if ((kTerm >= 1) && (kTerm <= 10))
YVEC = PSTC(kTerm, 1:nTM);
ZVEC = PSTC(kTerm, nTM+1:end);
alpha = 0.8;
plot3(XZ(1:nTM)-2, YVEC, ZVEC, ':', ...
'color', [0 0 0]+alpha);
plot3(XZ(1:nTM)+2, YVEC, ZVEC, ':', ...
'color', [0 0 0]+alpha);
endif;
YVEC = PSTC(kTerm+1, 1:iTM);
ZVEC = PSTC(kTerm+1, nTM+1:nTM+iTM);
plot3(XZ(1:iTM)-2, YVEC, ZVEC, 'k:');
plot3(XZ(1:iTM)+2, YVEC, ZVEC, 'k:'); % plot aggregated vectors (x=2 plane)
XVEC = [2 2 NaN 2 2];
for iPS=1:kTerm+1
iColor = mod(iPS-1,3) + 1; YVEC = [ ...
PSPC(iPS,iTM),
PSPC(iPS,iTM)-PSP(iPS,iTM),
NaN,
PSNC(iPS,iTM),
PSNC(iPS,iTM)-PSN(iPS,iTM)];
ZVEC = [ ...
PSPC(iPS,iTM+nTM),
PSPC(iPS,iTM+nTM)-PSP(iPS,iTM+nTM),
NaN,
PSNC(iPS,iTM+nTM),
PSNC(iPS,iTM+nTM)-PSN(iPS,iTM+nTM)];
plot3(XVEC, YVEC, ZVEC, CT(iColor), ...
'linewidth', 2);
endfor; % plot '+' vector addition cross-hair
iPS = kTerm + 1;
YVEC = [ ...
PSTC(iPS,iTM)-dCross/2,
PSTC(iPS,iTM)+dCross/2,
NaN,
PSTC(iPS,iTM),
PSTC(iPS,iTM)];
ZVEC = [ ...
PSTC(iPS,iTM+nTM),
PSTC(iPS,iTM+nTM),
NaN,
PSTC(iPS,iTM+nTM)-dCross/2,
PSTC(iPS,iTM+nTM)+dCross/2];
plot3(XVEC-4, YVEC, ZVEC, 'k-', ...
'linewidth', 1);
plot3(XVEC, YVEC, ZVEC, 'k-', ...
'linewidth', 1);
plot3( ...
XDOT+2, ...
YDOT+PSTC(iPS,iTM)/2, ...
ZDOT+PSTC(iPS,iTM+nTM)/2, ...
'k-', ...
'linewidth', 1); % print scene's frames
hold off;
fname = sprintf("./cfsanim3B/frame%05d.png", ...
iFrame);
print("-dpng", "-r400", "-Farial:10", fname);
iFrame++;
endfor; % iT=1:nT
endfor; % for kTerm=NTERM% plot _pceo.png (wher PCEO is Plot Coefficient Error in Octave)
figure(4);
X = I_CMAG4(:,1);
Y = diff(I_CMAG4(:,2:5), 1, 2); % diff of 1 in horizontal direction
plot(X, Y);
print("-dpng", "-r400", "-Farial:10", "_pceo.png");
I_CMAG4(nTermMax-4:nTermMax+6,:)% print 3 seconds of trailer still frames
for i=1:60
fname = sprintf("./cfsanim3B/frame%05d.png", ...
iFrame);
print("-dpng", "-r400", "-Farial:10", fname);
iFrame++;
endfor;% get final snapshot of side view
PV = view(); % Plot View
view([0 0]);
drawnow();
print("-dpng", "-r400", "-Farial:10", 'cfs3b3.png');
view(PV);% /tool/ffmpeg-4.3.1/bin/ffmpeg -f image2 -i frame%05d.png -vf scale=1080:-2,setsar=1:1 -f mp4 -q:v 0 -vcodec mpeg4 -r 20 _cfs3b.mp4
% /tool/ffmpeg-4.3.1/bin/ffmpeg -i _cfs3b.mp4 -c:v libvpx-vp9 -crf 30 -b:v 4M cfs3b.webm
% see https://trac.ffmpeg.org/wiki/Encode/VP9
Appendix B. Julia source code
The following Julia code generates the frames for the animation of the complex Fourier series tracing the letter ‘e’.
# file: cfse.jl (complex Fourier series tracing 'e')
# import Pkg; Pkg.add("Plots"); Pkg.add("PyPlot")
using Plots, Printf
pyplot()# cfse: Complex Fourier Series animation tracing the letter 'e'
# - media_mask bit 0: set to disable generation of frames as .png images
# - media_mask bit 1: set to disable generation of audio as .wav file(s)
# - iProc: 1-based index of CPU core
# - nProc: number of CPU cores
# example: cfse(1) # skips image frame generation
function cfse(media_mask::Int64=0, iProc=1, nProc=1)
# define some plotting parameters
nT = Int32(801) # number of Time samples
nTM = nT - 1 # number of Time sample Midpoints
T = range(0, stop=1, length=nT) # Time samples
nT2 = Int32(nTM/2)
nT4 = Int32(nTM/4)
dT = 1/nTM # difference between Time samples
TM = T[1:nTM] .+ dT/2 # Time sample Midpoints
THETAM = 2*pi*TM # angle samples
XZ = zeros(nT,1) # X Zeros
CSEQ = [:red, :green, :blue, :yellow, :white, :black] # Color SEQuence
trapScale = 0.5 # trapezoid scale
dCross = 0.15 # diameter of '+' cross-hair
yOffset = 0.5 # y offset of F(t)
nAudioBuf = 3530
nAng = 15
ANGDOT = range(0, stop=2*pi, length=nAng)
XDOT = zeros(nAng, 1)
YDOT = dCross/pi*cos.(ANGDOT)
ZDOT = dCross/pi*sin.(ANGDOT)
# parameterize function F for drawing letter 'e'
wChar = 0.1
rInner = 0.5
rOuter = rInner+wChar
cy = 0.4
thetaChar = 7*pi/4
F = zeros(1, 2*nT) # left half real, right half imaginary
M = zeros(2,1) # slope (Re and Im)
P = zeros(1,7) # Perimeter segments
P[1] = cy + rInner
P[2] = rInner * thetaChar
P[3] = wChar
P[4] = rOuter * thetaChar
P[5] = wChar
P[6] = cy + rOuter
P[7] = wChar
pSum = sum(P)
TS = P/pSum # Time Segments
TSC = cumsum(TS, dims=2) # Time Segments Cumulative
V = zeros(2, 7) # Vertices
V[1,1] = cy + rInner
V[2,1] = 0
V[1,2] = -rInner + rInner * cos.(thetaChar)
V[2,2] = rInner * sin.(thetaChar)
V[1,3] = wChar * cos.(thetaChar)
V[2,3] = wChar * sin.(thetaChar)
V[1,4] = -rOuter * cos.(thetaChar) + rOuter
V[2,4] = -rOuter * sin.(thetaChar)
V[1,5] = 0
V[2,5] = -wChar
V[1,6] = -cy - rOuter
V[2,6] = 0
V[1,7] = 0
V[2,7] = wChar
VC = cumsum(V, dims=2) # Vertices Cumulative
iSegment = 1 # segment 1
i0 = Int32(1)
tOffset = 0
nTS = 1 + Int32(floor(TS[1]/dT))
i1 = i0 + nTS - 1
M = V[:,1] ./ TS[1]
F[i0:i1] = M[1] * T[i0:i1]
F[nT+i0:nT+i1] .= 0
iSegment += 1 # segment 2
i0 = i1 + 1
tOffset = T[i0] - TSC[1]
nTS = 1 + Int32(floor((TS[2] - tOffset)/dT))
i1 = i0 + nTS - 1
period = 8*TS[2]/7
THETA = 2*pi/period*(T[i0:i1] .- TSC[1])
F[i0:i1] = cy .+ rInner*cos.(THETA)
F[nT+i0:nT+i1] = rInner*sin.(THETA)
iSegment += 1 # segment 3
i0 = i1 + 1
tOffset = T[i0] - TSC[2]
nTS = 1 + Int32(floor((TS[3] - tOffset)/dT))
i1 = i0 + nTS - 1
M = V[:,3] ./ TS[3]
F[i0:i1] = VC[1, 2] .+ M[1] .* (T[i0:i1] .- TSC[2])
F[nT+i0:nT+i1] = VC[2, 2] .+ M[2] * (T[i0:i1] .- TSC[2])
iSegment += 1 # segment 4
i0 = i1 + 1
tOffset = T[i0] - TSC[3]
nTS = 1 + Int32(floor((TS[4] - tOffset)/dT))
i1 = i0 + nTS - 1
period = 8*TS[4]/7
THETA = thetaChar .- 2*pi/period*(T[i0:i1] .- TSC[3])
F[i0:i1] = cy .+ rOuter*cos.(THETA)
F[nT+i0:nT+i1] = rOuter*sin.(THETA)
iSegment += 1 # segment 5
i0 = i1 + 1
tOffset = T[i0] - TSC[4]
nTS = 1 + Int32(floor((TS[5] - tOffset)/dT))
i1 = i0 + nTS - 1
M = V[:,5] ./ TS[5]
F[i0:i1] .= VC[1, 4]
F[nT+i0:nT+i1] = VC[2, 4] .+ M[2] * (T[i0:i1] .- TSC[4])
iSegment += 1 # segment 6
i0 = i1 + 1
tOffset = T[i0] - TSC[5]
nTS = 1 + Int32(floor((TS[6] - tOffset)/dT))
i1 = i0 + nTS - 1
M = V[:,6] ./ TS[6]
F[i0:i1] = VC[1, 5] .+ M[1] * (T[i0:i1] .- TSC[5])
F[nT+i0:nT+i1] .= VC[2, 5]
iSegment += 1 # segment 7
i0 = i1 + 1
tOffset = T[i0] - TSC[6]
nTS = 1 + Int32(floor((TS[7] - tOffset)/dT))
i1 = i0 + nTS - 1
M = V[:,7] ./ TS[7]
F[i0:i1] .= VC[1, 6]
F[nT+i0:nT+i1] = VC[2, 6] .+ M[2] * (T[i0:i1] .- TSC[6])
# apply y offset shift then take midpoints
F[1:nT] .+= yOffset
FM = [(F[1:nT-1]+F[2:nT])/2 (F[nT+1:end-1]+F[nT+2:end])/2]
# integrate complex coefficient
NTERM = [1:10; 12:2:22]
nTerm = size(NTERM,1)
nTermMax = last(NTERM)
FEM = zeros(1, 2*nTM);
CR = zeros(5, 2*nTermMax+1)
CI = zeros(5, 2*nTermMax+1)
for iTerm=-nTermMax:nTermMax
iA = iTerm + nTermMax + 1
EM = [cos.(-iTerm*THETAM) sin.(-iTerm*THETAM)]
FEM[1:nTM] = FM[1:nTM] .* EM[1:nTM] -
FM[nTM+1:end] .* EM[nTM+1:end]
FEM[nTM+1:end] = FM[1:nTM] .* EM[nTM+1:end] +
FM[nTM+1:end] .* EM[1:nTM]
V = reshape(FEM .* dT, (nTM, 2))
VR = reshape(V[:,1], (nT4, 4))
CR[1:4,iA] = sum(VR, dims=1)'
VI = reshape(V[:,2], (nT4, 4))
CI[1:4,iA] = sum(VI, dims=1)'
end
CR[5,:] = sum(CR[1:4,:], dims=1)
CI[5,:] = sum(CI[1:4,:], dims=1)
if (iProc == 1)
fp = open("ecoeff.dat", "w")
write(fp, CR[5,:])
write(fp, CI[5,:])
close(fp)
end
# calculate common ratio r that best fits coefficients
A = zeros(2, 1) # (a ccw terms; a cw terms)
R = zeros(2, 1) # [r ccw terms, r cw terms]
I_CMAG4 = zeros(2*nTermMax+1, 5) # columns 3 thru 5 are for testing
I_CMAG4[:,1] .= -nTermMax:nTermMax
I_CMAG4[:,2] .= (CR[5,:].^2 + CI[5,:].^2).^0.5
X = [ones(nTermMax,1) I_CMAG4[nTermMax+2:end,1]]
Y = log.(I_CMAG4[nTermMax+2:end,2])
PHI = (X'*X) \ (X'*Y) # is faster than PHI = inv(X'*X)*X'*Y
EPHI = exp.(PHI)
A[1] = EPHI[1]
R[1] = EPHI[2]
X = [ones(nTermMax,1) abs.(I_CMAG4[1:nTermMax,1])]
Y = log.(I_CMAG4[1:nTermMax,2])
PHI = (X'*X) \ (X'*Y) # is faster than PHI = inv(X'*X)*X'*Y
EPHI = exp.(PHI)
A[2] = EPHI[1]
R[2] = EPHI[2]
if (nProc == 1)
display("R"); display(R)
end
# plot L.S.E. fit to natural log of coefficient magnitudes
plot(I_CMAG4[nTermMax+2:end,1],
log.(A[1]*R[1].^I_CMAG4[nTermMax+2:end,1]),
color = :black);
plot!(I_CMAG4[1:nTermMax,1],
log.(A[2]*R[2].^abs.(I_CMAG4[1:nTermMax,1])),
color = :black);
ICOLOR = mod.(abs.(I_CMAG4[:,1]),3) .== 0;
plot!(I_CMAG4[ICOLOR,1],
log.(I_CMAG4[ICOLOR,2]),
seriestype = :scatter,
markersize = 8,
color = :red);
ICOLOR = mod.(abs.(I_CMAG4[:,1]),3) .== 1;
plot!(I_CMAG4[ICOLOR,1],
log.(I_CMAG4[ICOLOR,2]),
seriestype = :scatter,
markersize = 8,
color = :green);
ICOLOR = mod.(abs.(I_CMAG4[:,1]),3) .== 2;
plot!(I_CMAG4[ICOLOR,1],
log.(I_CMAG4[ICOLOR,2]),
seriestype = :scatter,
markersize = 8,
color = :blue);
strTitleCoefficient = @sprintf(
"%s\n+ index common ratio r = %.4f\n- index common ratio r = %.4f",
"LSE fit to coefficient magnitude data",
R[1],
R[2]);
plot!(legend = false,
xlabel = "coefficient index n",
ylabel = "ln(coefficient magnitude)",
title = strTitleCoefficient);
savefig("_pc.png")
# plot coefficient magnitudes as trapezoids
SMB = zeros(2, 1) # Sum M*B; dims 1: +/- coefficients
SB = zeros(2, 1) # Sum B; dims 1: +/- coefficients
TH = zeros(2, nTermMax+1) # Trapezoid Height; dims 1: +/- coefficients; dims2: trapezoid height (1-based coefficient index)
X = zeros(2, 2) # dims 1: +/- coefficients; dims 2: x of tall/short height
M = zeros(2, 1) # trapezoid slope; dims 1: +/- coefficients
B = zeros(2, 1) # trapezoid base width; dims 1: +/- coefficients
X[:,2] = 2 .* (1 .- R .^ 0.5)
TH[:,1] = [I_CMAG4[nTermMax+1,2], I_CMAG4[nTermMax+1,2]] ./ 2 # split DC component into +/- coefficients
TH[:,1] = TH[:,1] ./ (1 .- R) # trapezoid height: dims 1: +/- coefficients
XTRAP = [
X[1,1], X[1,1], X[1,2], X[1,2], X[1,1],
NaN,
-X[2,1], -X[2,1], -X[2,2], -X[2,2], -X[2,1]]
YTRAP = [
0, TH[1,1], TH[1,1]*R[1]^0.5, 0, 0,
NaN,
0, TH[2,1], TH[2,1]*R[2]^0.5, 0, 0]
plot(XTRAP, YTRAP,
color = CSEQ[1], # coefficient 0
legend = false,
xlims = (-2,2),
xticks = -2:2:2)
# SB += B # intentionally skip inclusion of DC component
# SMB += M .* B # in calculation of weighted average slope
I_CMAG4[nTermMax+1,3] = sum(TH[:,1] .* (1 .- R), dims=1)[1]
for jTerm=1:nTermMax
iColor::Int16 = mod(jTerm,3) + 1;
X[:,1] = X[:,2]
X[:,2] = 2 .* (1 .- R .^ ((jTerm+1)/2))
B = abs.(diff(X, dims=2))
TH[:,jTerm+1] = [
I_CMAG4[nTermMax+1+jTerm,2];
I_CMAG4[nTermMax+1-jTerm,2]] ./
(R .^ (jTerm/2) .* (1 .- R))
M = TH[:,jTerm+1] .* (R .- 1) ./ B
XTRAP = [
X[1,1], X[1,1], X[1,2], X[1,2], X[1,1],
NaN,
-X[2,1], -X[2,1], -X[2,2], -X[2,2], -X[2,1]]
YTRAP = [
0, TH[1,jTerm+1], TH[1,jTerm+1]*R[1]^0.5, 0, 0,
NaN,
0, TH[2,jTerm+1], TH[2,jTerm+1]*R[2]^0.5, 0, 0]
plot!(XTRAP, YTRAP,
color = CSEQ[iColor]) # coefficient +/- jTerm
SB += B
SMB += M .* B
I_CMAG4[nTermMax+1+jTerm,3] = TH[1,jTerm+1] * R[1]^(jTerm/2) * (1 - R[1])
I_CMAG4[nTermMax+1-jTerm,3] = TH[2,jTerm+1] * R[2]^(jTerm/2) * (1 - R[2])
end
MWA = SMB ./ SB # M (i.e., slope) Weighted Average
plot!(
[0, 2, NaN, 0, -2],
[-MWA[1]*2, 0, NaN, -MWA[2]*2, 0],
linestyle = :dot,
color = :orange)
savefig("_pt.png")
if (nProc == 1)
display("MWA"); display(MWA)
end
# calculate partial series terms
PS = zeros(2*nTM, 2*nTermMax+1) # Partial Sum
iTerm0 = nTermMax + 1
PS[1:nTM,1] .= CR[5,iTerm0]
PS[nTM+1:end,1] .= CI[5,iTerm0]
for jTerm=1:nTermMax
iPS = 2 * jTerm
iTerm = iTerm0 + jTerm
EM = [cos.(jTerm*THETAM) sin.(jTerm*THETAM)] # nTM x 2
PS[1:nTM,iPS] =
CR[5,iTerm] .* EM[:,1] -
CI[5,iTerm] .* EM[:,2]
PS[nTM+1:end,iPS] =
CI[5,iTerm] .* EM[:,1] +
CR[5,iTerm] .* EM[:,2]iPS = 2 * jTerm + 1
iTerm = iTerm0 - jTerm
EM = [cos.(jTerm*THETAM) sin.(-jTerm*THETAM)] # nTM x 2
PS[1:nTM,iPS] =
CR[5,iTerm] .* EM[:,1] -
CI[5,iTerm] .* EM[:,2]
PS[nTM+1:end,iPS] =
CI[5,iTerm] .* EM[:,1] +
CR[5,iTerm] .* EM[:,2]
end
PSC = cumsum(PS, dims=2);
# split terms according to direction of rotation
PSP = zeros(2*nTM, nTermMax+1) # Partial Sum Positive rotation
PSP[:,1] = PS[:,1] ./ 2
PSP[:,2:end] = PS[:,2:2:2*nTermMax+1]
PSPC = cumsum(PSP, dims=2) # Partial Sum Positive rotation Cumulative
PSN = zeros(2*nTM, nTermMax+1) # Partial Sum Negative rotation
PSN[:,1] = PS[:,1] ./ 2
PSN[:,2:end] = PS[:,3:2:2*nTermMax+1]
PSNC = cumsum(PSN, dims=2) # Partial Sum Negative rotation Cumulative
PST = PSP + PSN # Partial Sum Total
PSTC = cumsum(PST, dims=2) # Partial Sum Total Cumulative
# test PSP and PSN in column 4 of I_CMAG4
kTM = 100 # test should work for any value of kTM where 1<kTM<=nTM
I_CMAG4[nTermMax+1,4] = (PSP[kTM,1]^2 + PSP[kTM+nTM,1]^2)^0.5
I_CMAG4[nTermMax+1,4] += (PSN[kTM,1]^2 + PSN[kTM+nTM,1]^2)^0.5
for jTerm=2:nTermMax+1
I_CMAG4[nTermMax+jTerm,4] = (PSP[kTM,jTerm]^2 + PSP[kTM+nTM,jTerm]^2)^0.5
I_CMAG4[nTermMax+2-jTerm,4] = (PSN[kTM,jTerm]^2 + PSN[kTM+nTM,jTerm]^2)^0.5
end
# calculate frame offsets where new terms are added to partial series
FRAMEOFFSET = nTM/2 * ones(Int32, nTerm+1, 1)
FRAMEOFFSET[1] = 0
FRAMEOFFSET[nTerm+1] *= 2 # n=24 is slow motion
FRAMEOFFSET = cumsum(FRAMEOFFSET, dims=1)
# generate audio tracks
# https://isip.piconepress.com/projects/speech/software/tutorials/production/fundamentals/v1.0/section_02/s02_01_p05.html
if ((media_mask&0x2==0) && (iProc==1))
io = open("./ez2armst/_ez2armst.wav", "w")
# write 44 byte .WAV header
write(io, "RIFF") # RIFF tag (4 bytes)
uival::UInt32 = 0 # (to be completed after writing of audio data)
write(io, uival) + # tag block size (4 bytes)
write(io, "WAVE") + # WAVE tag (4 bytes)
write(io, "fmt ") # format description header (4 bytes)
uival = 16
write(io, uival) # tag block size (4 bytes)
usval::UInt16 = 1
write(io, usval) # specify PCM encoding (2 bytes)
usval = 2
write(io, usval) # specify stereo sampling (2 bytes)
uival = 44100
write(io, uival) # specify sampling rate (4 bytes)
uival *= 4
write(io, uival) # specify bytes/second (4 bytes)
usval = 4
write(io, usval) # specify block alignment (2 bytes)
usval = 16
write(io, usval) # specify bits per sample (2 bytes)
write(io, "data") # data description header (4 bytes)
uival = 0 # (to be completed after writing of audio data)
write(io, uival) # specify # bytes of audio data (4 bytes)
# write audio track data (subtracting off DC offset in PSC[1,:]
AUDIO = zeros(Int16, nT4, 2)
AUDIO2 = zeros(Int16, nT4, 2)
maxPSC = maximum(abs.(PSC - repeat(PSC[:,1], outer=(1,2*nTermMax+1))))
if (nProc == 1)
display(maxPSC)
end
I = 1:4:nTM
J = 1:2:nT2
for iTerm = 1:nTerm
kTerm = NTERM[iTerm]
iPS = 2*kTerm + 1 # +1 to skip DC component in PS
if (iTerm != nTerm)
AUDIO[:,1] = Int16.(round.((PSC[I,iPS]-PSC[I,1]) .*
(0x7fff*0.8/maxPSC)))
AUDIO[:,2] = Int16.(round.((PSC[I.+nTM,iPS]-PSC[I.+nTM,1]) .*
(0x7fff*0.8/maxPSC)))
for iBuf=1:nAudioBuf
write(io, AUDIO'[:])
end
else
AUDIO[:,1] = Int16.(round.((PSC[J,iPS]-PSC[J,1]) .*
(0x7fff*0.8/maxPSC)))
AUDIO[:,2] = Int16.(round.((PSC[J.+nTM,iPS]-PSC[J.+nTM,1]) .*
(0x7fff*0.8/maxPSC)))
AUDIO2[:,1] = Int16.(round.((PSC[J.+nT2,iPS]-PSC[J.+nT2,1]) .*
(0x7fff*0.8/maxPSC)))
AUDIO2[:,2] = Int16.(round.((PSC[J.+(nTM+nT2),iPS]-PSC[J.+(nTM+nT2),1]) .*
(0x7fff*0.8/maxPSC)))
for iBuf=1:nAudioBuf
write(io, AUDIO'[:])
write(io, AUDIO2'[:])
end
end
end # for iTerm=1:nTerm
# complete .WAV file header and close file
fsize::UInt32 = position(io)
seek(io, 4)
write(io, fsize - UInt32(8))
seek(io, 40)
write(io, fsize - UInt32(44))
close(io)
end
# animate
if (media_mask&0x1 == 0)
if (nProc == 1)
print("generating 3D animation frames ") # start progress report line
end
iTerm1 = 1 + div((iProc-1)*nTerm, nProc, RoundNearest)
iTerm2 = div(iProc*nTerm, nProc, RoundNearest)
if (nProc > 1)
str = @sprintf("partition %d:%d ", iTerm1, iTerm2)
print(str)
end
for iTerm = iTerm1:iTerm2 # partitioned loop
# for iTerm = 1:nTerm # original unpartitioned loop
print(".") # progress report: one '.' for each processed term in NTERM
kTerm = NTERM[iTerm]
iFrame::Int32 = FRAMEOFFSET[iTerm]
for iTM = 1:nTM
# typically skip every other frame for speed
if (iTerm != nTerm)
if (mod(iTM,2)==0)
continue;
end
end
# restore background plot but with new title
X = 2*ones(nT, 1);
plot(X,F[1:nT],F[nT+1:2*nT],
color = :cyan)
plot!( # M (i.e., slope) Weighted Average
[0, 2, NaN, 0, -2],
[-trapScale*MWA[1]*2, 0, NaN, -trapScale*MWA[2]*2, 0],
[0, 0, NaN, 0, 0],
color = :orange,
linestyle = :dot,
linewidth = 0.5);
strTitle = @sprintf("%s\n(%s=%.2f, %s=%4d,\nt=%.4f, n=+/-%d)",
"complex Fourier series tracing letter 'e'",
"trapezoid height scale",
trapScale,
"sound frequency scale",
nAudioBuf,
iTM/nTM,
kTerm)
plot!(background=:white,
camera = (73, 20),
xlabel = "trapezoid baseline ",
ylabel = "Re Sn",
zlabel = "\n\nIm Sn",
title = strTitle);
plot!(legend = false,
xlims = (-2,2),
ylims = (-0.99,1.99),
zlims = (-1,1),
xticks = 0:2:2,
yticks = -1:1:2,
zticks = -1:1:1,
figsize = (800,450))
plot!([-2, 2], [0, 0], [0, 0],
color = RGB(0.82,0.82,0.82))
plot!(-X,F[1:nT],F[nT+1:2*nT], color="cyan")
# plot zigzag aggregated vectors (x=-2 plane)
XVEC = [-2, -2]
YVEC = [0, PS[iTM,1]]
ZVEC = [0, PS[iTM+nTM,1]]
plot!(XVEC, YVEC, ZVEC,
color = CSEQ[1],
linewidth = 1);
XVEC = [-2, -2, -2]
for iPS=2:2:2*kTerm
iColor::Int16 = mod(floor(iPS/2),3) + 1;
YVEC = [PSC[iTM,iPS-1], PSC[iTM,iPS], PSC[iTM,iPS+1]]
ZVEC = [PSC[iTM+nTM,iPS-1], PSC[iTM+nTM,iPS], PSC[iTM+nTM,iPS+1]]
plot!(XVEC, YVEC, ZVEC,
color = CSEQ[iColor],
linewidth = 1);
end
# plot axes
plot!(xlims(), [0, 0], [0, 0],
color = :gray85);
plot!([2, 2], [0, ylims()[2]], [0, 0],
color = :gray85);
# plot trapezoids
for iPS=1:kTerm+1
iColor::Int16 = mod(iPS-1,3) + 1;
S = R .^ ((iPS-1)/2) .* (1 .- R)
# plot trapezoids (mid planes)
XTRAP = [
2*(1-R[1]^((iPS-1)/2)),
2*(1-R[1]^((iPS-1)/2)),
2*(1-R[1]^(iPS/2)),
2*(1-R[1]^(iPS/2)),
2*(1-R[1]^((iPS-1)/2)),
NaN,
-2*(1-R[2]^((iPS-1)/2)),
-2*(1-R[2]^((iPS-1)/2)),
-2*(1-R[2]^(iPS/2)),
-2*(1-R[2]^(iPS/2)),
-2*(1-R[2]^((iPS-1)/2))]
YTRAP = [
0,
trapScale*PSP[iTM,iPS]/S[1],
trapScale*PSP[iTM,iPS]/S[1] * R[1]^0.5,
0,
0,
NaN,
0,
trapScale*PSN[iTM,iPS]/S[2],
trapScale*PSN[iTM,iPS]/S[2] * R[2]^0.5,
0,
0]
ZTRAP = [
0,
trapScale*PSP[iTM+nTM,iPS]/S[1],
trapScale*PSP[iTM+nTM,iPS]/S[1] * R[1]^0.5,
0,
0,
NaN,
0,
trapScale*PSN[iTM+nTM,iPS]/S[2],
trapScale*PSN[iTM+nTM,iPS]/S[2] * R[2]^0.5,
0,
0]
plot!(XTRAP, YTRAP, ZTRAP,
color = CSEQ[iColor],
linewidth = 0.3);
## fill test column 5 of I_CMAG4 ##
# left edge of trapezoid is
# x(i) = 2*(1 - r^(i/2))
#
# width of trapezoid base is
# b(i) = x(i+1) - x(i)
# = 2*r^(i/2)*(1 - r^(1/2))
#
# given lefth height (h(i)), right height of trapezoid is
# h(i+1) = h(i)*r^(1/2)
#
# trapezoid area is
# |C| = h(i)*b(i) - 1/2(h(i) - h(i+1))*b(i)
# = h(i)*r^(i/2)*(1 - r)
#
# where h is trapezoid height, r is common ratio,
# and i is the power of the term represented by
# the trapezoid area
if ((kTerm == nTermMax) && (iTM == 1))
I = [2; 8]
H = ((YTRAP[I].^2 + ZTRAP[I].^2).^0.5) ./ trapScale
if (iPS > 1)
I = [nTermMax+iPS; nTermMax+2-iPS]
I_CMAG4[I,5] = H .* R .^ ((iPS-1)/2) .* (1 .- R)
else
ADC = H .* (1 .- R) # Area of DC trapezoids
I_CMAG4[nTermMax+1,5] = sum(ADC, dims=1)[1]
end
end
end
# plot dotted trail of partial sum
if ((kTerm >= 1) && (kTerm <= 10))
YVEC = PSTC[1:nTM,kTerm]
ZVEC = PSTC[nTM+1:end,kTerm]
plot!(XZ[1:nTM] .- 2, YVEC, ZVEC,
linestyle = :dot,
color = :gray80)
plot!(XZ[1:nTM] .+ 2, YVEC, ZVEC,
linestyle = :dot,
color = :gray80)
end
YVEC = PSTC[1:iTM,kTerm+1]
ZVEC = PSTC[nTM+1:nTM+iTM,kTerm+1]
plot!(XZ[1:iTM] .- 2, YVEC, ZVEC,
linestyle = :dot,
color = :black)
plot!(XZ[1:iTM] .+ 2, YVEC, ZVEC,
linestyle = :dot,
color = :black)
# plot 2-arm aggregated vectors (x=2 plane)
XVEC = [2, 2, NaN, 2, 2]
for iPS=1:kTerm+1
iColor::Int16 = mod(iPS-1,3) + 1
YVEC = [PSPC[iTM,iPS],
PSPC[iTM,iPS]-PSP[iTM,iPS],
NaN,
PSNC[iTM,iPS],
PSNC[iTM,iPS]-PSN[iTM,iPS]]
ZVEC = [PSPC[iTM+nTM,iPS],
PSPC[iTM+nTM,iPS]-PSP[iTM+nTM,iPS],
NaN,
PSNC[iTM+nTM,iPS],
PSNC[iTM+nTM,iPS]-PSN[iTM+nTM,iPS]]
plot!(XVEC, YVEC, ZVEC,
color = CSEQ[iColor],
linewidth = 2)
end
# plot '+' vector addition cross-hair
iPS = kTerm + 1
YVEC = [
PSTC[iTM,iPS]-dCross/2,
PSTC[iTM,iPS]+dCross/2,
NaN,
PSTC[iTM,iPS],
PSTC[iTM,iPS]]
ZVEC = [
PSTC[iTM+nTM,iPS],
PSTC[iTM+nTM,iPS],
NaN,
PSTC[iTM+nTM,iPS]-dCross/2,
PSTC[iTM+nTM,iPS]+dCross/2];
plot!(XVEC .- 4, YVEC, ZVEC,
color = :black,
linewidth = 1);
plot!(XVEC, YVEC, ZVEC,
color = :black,
linewidth = 1);
plot!(XDOT .+ 2,
YDOT .+ PSTC[iTM,iPS]/2,
ZDOT .+ PSTC[iTM+nTM,iPS]/2,
color = :black,
linewidth = 1);
# print scene's frames
iFrame += 1
fname = @sprintf("./ez2armst/frame%05d.png", iFrame);
savefig(fname)
# print("-dpng", "-r400", "-Farial:10", fname);
end # iT=1:nT
end # for kTerm=NTERM
print("\n") # end progress report line
end # if (media_mask&0x1 == 0)
# plot _pce.png (where PCE is Plot Coefficient Error) to show
# consistency of calculated coefficient magnitude |C| values in I_CMAG4
if (nProc == 1)
print("plotting (_pce.png) the numerical (e.g., trig approximation) errors in coefficient calculations ...\n")
X = I_CMAG4[:,1]
Y = diff(I_CMAG4[:,2:5], dims=2)
plot(X, Y)
savefig("_pce.png")
# display("I_CMAG4[nTermMax+1-5:nTermMax+1+5,:]"); display(I_CMAG4[nTermMax+1-5:nTermMax+1+5,:])
end
# /tool/ffmpeg-4.3.1/bin/ffmpeg -f image2 -i frame%05d.png -i _ez2armst.mp3 -vf scale=1080:-2,setsar=1:1 -f mp4 -q:v 0 -vcodec mpeg4 -r 20 _ez2armst.mp4
# /tool/ffmpeg-4.3.1/bin/ffmpeg -f image2 -i frame%05d.png -vf scale=1080:-2,setsar=1:1 -f mp4 -q:v 0 -vcodec mpeg4 -r 20 _ez2armst.mp4
# /tool/ffmpeg-4.3.1/bin/ffmpeg -i _ez2armst.mp4 -c:v libvpx-vp9 -crf 30 -b:v 4M _ez2armst.webm
# see https://trac.ffmpeg.org/wiki/Encode/VP9
end
Appendix C. Julia source code
The following Julia code generates the animation frames of the complex geometric series.
# import Pkg; Pkg.add("Plots"); Pkg.add("PyPlot")
using Plots, Printf
pyplot()# cgs: Complex Geometric Series
# isTrapezoid: 0/1 flag for rotating triangles/trapezoids
# a: coefficient
# r: common ratio
# n: power of last term in partial series
# nAudioBuf: audio buffer count per animated fundamental period
# example: cgs()
function cgs(isTrapezoid::Int64=0,
a::Float64=1.0,
r::Float64=0.5,
n::Int64=10,
nAudioBuf::Int64=3530)
# define some plotting parameters
nT = Int32(400) # # of non-repeated Time samples
nT1 = nT + 1 # # of Time samples
T = range(0, stop=1, length=nT1) # Time samples
nT2 = Int32(nT/2)
THETA = 2*pi*T # angle samples of plot's 'o' symbol
NTERM = 1:n
nTerm = size(NTERM,1)
nTermMax = last(NTERM)
maxSum = a/(1-abs(r)) # for setting plot limits
maxRadius = a*abs(r)/(1-r^2)
dCross = 0.15 # diameter of '+' symbol
trapScale = 0.5 # trapezoid scale
CSEQ = [:red, :green, :blue, :yellow, :white, :black] # Color SEQuence
# calculate partial series terms
PS = zeros(2*nT, nTermMax+1) # Partial Sum
PS[1:nT,1] .= a
for jTerm=1:nTermMax
EM = [ sin.(jTerm*THETA)] # nT x 2
PS[1:nT,jTerm+1] = (a * r^jTerm) .*
cos.(jTerm*THETA[1:end-1])
PS[nT+1:end,jTerm+1] = (a * r^jTerm) .*
sin.(jTerm*THETA[1:end-1])
end
PSC = cumsum(PS, dims=2);
# calculate frame offsets where new terms are added to partial series
FRAMEOFFSET = nT * ones(Int32, nTerm+1, 1)
FRAMEOFFSET[1] = 0
FRAMEOFFSET = cumsum(FRAMEOFFSET, dims=1)
# generate audio tracks
io = open("./cgs/cgs.wav", "w")
# write 44 byte .WAV header
write(io, "RIFF") # RIFF tag (4 bytes)
uival::UInt32 = 0 # (to be completed after writing of audio data)
write(io, uival) + # tag block size (4 bytes)
write(io, "WAVE") + # WAVE tag (4 bytes)
write(io, "fmt ") # format description header (4 bytes)
uival = 16
write(io, uival) # tag block size (4 bytes)
usval::UInt16 = 1
write(io, usval) # specify PCM encoding (2 bytes)
usval = 2
write(io, usval) # specify stereo sampling (2 bytes)
uival = 44100
write(io, uival) # specify sampling rate (4 bytes)
uival *= 4
write(io, uival) # specify bytes/second (4 bytes)
usval = 4
write(io, usval) # specify block alignment (2 bytes)
usval = 16
write(io, usval) # specify bits per sample (2 bytes)
write(io, "data") # data description header (4 bytes)
uival = 0 # (to be completed after writing of audio data)
write(io, uival) # specify # bytes of audio data (4 bytes)
# write audio track data (subtracting off DC offset in PSC[1,:]
AUDIO = zeros(Int16, nT2, 2)
I = 1:2:nT
for iTerm = 1:nTerm
kTerm = NTERM[iTerm]
if (mod(iTerm,2) == 0) # if power of last term is even
minmaxmid = a*(1 - (abs(r))^(kTerm+2)) / (1 - (abs(r))^2)
else # else power of last term is odd
minmaxmid = a*(1 - (abs(r))^(kTerm+1)) / (1 - (abs(r))^2)
end
AUDIO[:,1] = Int16.(round.((PSC[I,kTerm+1] .- minmaxmid) .*
(0x7fff*0.8/maxRadius)))
AUDIO[:,2] = Int16.(round.(PSC[I.+nT,kTerm+1] .*
(0x7fff*0.8/maxRadius)))
for iBuf=1:nAudioBuf
write(io, AUDIO'[:])
end
end # for iTerm=1:nTerm
# complete .WAV file header and close file
fsize::UInt32 = position(io)
seek(io, 4)
write(io, fsize - UInt32(8))
seek(io, 40)
write(io, fsize - UInt32(44))
close(io)
# define circle traced by complex geometric series
ycenter = a / (1 - r^2)
zcenter = 0
TH = range(0, stop=2*pi, length=nT1)
XCIRC = 2*ones(nT1, 1)
YCIRC = ycenter .+ maxRadius .* cos.(TH)
ZCIRC = zcenter .+ maxRadius .* sin.(TH)
# animate
print("generating 3D animation frames ") # start progress report line
PO = Vector{Plots.Plot{Plots.PyPlotBackend}}(undef,nTerm) # Plot Objects
for iTerm = 1:nTerm
print(".") # progress report: one '.' for each processed term in NTERM
kTerm = NTERM[iTerm]
iFrame::Int32 = FRAMEOFFSET[iTerm]
for iT = 1:nT
# restore background plot but with new title
if (isTrapezoid > 0)
strTitle = @sprintf(
"%s\n(%s=%.2f, %s=%.2fexp(iθ),\n%s=%.2f, %s=%d,\nt=%.4f, n=%d)",
"complex geometric series",
"coefficient a", a,
"common ratio r", r,
"trapezoid scale", trapScale,
"sound frequency scale", nAudioBuf,
(iT-1)/nT,
kTerm)
else
strTitle = @sprintf(
"%s\n(%s=%.2f, %s=%.2fexp(iθ),\n%s=%d,\nt=%.4f, n=%d)",
"complex geometric series",
"coefficient a", a,
"common ratio r", r,
"sound frequency scale", nAudioBuf,
(iT-1)/nT,
kTerm)
end
PO[iTerm] = plot(background=:white,
camera = (73, 20),
xlabel = "x",
ylabel = "y",
zlabel = "z",
title = strTitle)
plot!(legend = false,
xlims = (0,2),
ylims = (-0.5,2.0),
zlims = (-0.85,0.85),
xticks = 0:2:2,
yticks = 0:1:2,
zticks = -1:1:1,
figsize = (800,450))
# plot axes
plot!([2, 2], [0, ylims()[2]], [0, 0],
color = :gray85);
# plot trapezoids or triangles
if (isTrapezoid > 0)
# plot trapezoids
for iPS=1:kTerm+1
s = r ^ ((iPS-1)/2) * (1 - r)
iColor::Int16 = mod(iPS-1,3) + 1;
XTRAP = [
2*r^((iPS-1)/2),
2*r^((iPS-1)/2),
2*r^(iPS/2),
2*r^(iPS/2),
2*r^((iPS-1)/2)]
YTRAP = [
0,
trapScale*PS[iT,iPS]/s,
trapScale*PS[iT,iPS]/s * r^0.5,
0,
0]
ZTRAP = [
0,
trapScale*PS[iT+nT,iPS]/s,
trapScale*PS[iT+nT,iPS]/s * r^0.5,
0,
0]
plot!(XTRAP, YTRAP, ZTRAP,
color = CSEQ[iColor],
linewidth = 0.3);
end
else
# plot triangles
for iPS=1:kTerm+1
iColor::Int16 = mod(iPS-1,3) + 1;
XTRI = [
0,
2*r^((iPS-1)/2),
2*r^((iPS-1)/2),
0]
YTRI = [
0,
0,
2*PS[iT,iPS]/XTRI[3],
0]
ZTRI = [
0,
0,
2*PS[iT+nT,iPS]/XTRI[3],
0]
plot!(XTRI, YTRI, ZTRI,
color = CSEQ[iColor],
linewidth = 0.3);
end
end
plot!(XCIRC, YCIRC, ZCIRC, color="cyan")
# plot dotted trail of partial sum
if ((kTerm >= 1) && (kTerm <= 10))
YVEC = PSC[1:nT,kTerm]
ZVEC = PSC[nT+1:end,kTerm]
plot!(XCIRC[1:nT], YVEC, ZVEC,
linestyle = :dot,
color = :gray80)
end
YVEC = PSC[1:iT,kTerm+1]
ZVEC = PSC[nT+1:nT+iT,kTerm+1]
plot!(XCIRC[1:iT], YVEC, ZVEC,
linestyle = :dot,
color = :black)
# plot 1-arm aggregated vectors (x=2 plane)
XVEC = [2, 2]
for iPS=1:kTerm+1
iColor::Int16 = mod(iPS-1,3) + 1
YVEC = [
PSC[iT,iPS],
PSC[iT,iPS]-PS[iT,iPS]]
ZVEC = [
PSC[iT+nT,iPS],
PSC[iT+nT,iPS]-PS[iT+nT,iPS]]
plot!(XVEC, YVEC, ZVEC,
color = CSEQ[iColor],
linewidth = 2)
end
# print scene's frames
iFrame += 1
fname = @sprintf("./cgs/frame%05d.png", iFrame);
savefig(fname)
end # iT=1:nT
end # for kTerm=NTERM
print("\n") # end progress report line
# /tool/ffmpeg-4.3.1/bin/ffmpeg -f image2 -i frame%05d.png -i cgs.mp3 -vf scale=1080:-2,setsar=1:1 -f mp4 -q:v 0 -vcodec mpeg4 -r 20 cgs.mp4
# /tool/ffmpeg-4.3.1/bin/ffmpeg -f image2 -i frame%05d.png -vf scale=1080:-2,setsar=1:1 -f mp4 -q:v 0 -vcodec mpeg4 -r 20 cgs.mp4
# /tool/ffmpeg-4.3.1/bin/ffmpeg -i cgs.mp4 -c:v libvpx-vp9 -crf 30 -b:v 4M cgs.webm
# see https://trac.ffmpeg.org/wiki/Encode/VP9
end
Appendix D. Julia source code
The following Julia code distributes the complex Fourier series animation frame generation to six CPU cores and completed in 11.1 minutes which is a factor of 10 faster than the 120 minutes that GNU Octave took to generate the animation frames.
# file: cfsed.jl (complex Fourier series tracing 'e' distributed)
# import Pkg.add("Distributed")
using Distributedaddprocs(4, exeflags="--project=.")
@everywhere begin
include("cfse.jl")
endnProc = 8
pmap(1:nProc) do i
cfse(0, i, nProc)
end
Appendix E. Julia source code
The following is Julia code that generates miscellaneous plots including the geometric series self-similarity animation and a plot of the geometric series as the origin to the reference frame of the body of knowledge of mathematics.
# file: some1.jl (miscellaneous code for Summer Of Math Exposition #1)
# import Pkg; Pkg.add("Plots"); Pkg.add("PyPlot")
using Measures, Printf, LaTeXStrings, Plots
pyplot(dpi=600)function tstplt6()
a = 1:3
b = [1, 2, 7]
plot(a, b,
label = "randData",
xlabel = "numbers",
ylabel = "Rand data",
color = :red,
legend = :topleft,
tick_direction = :in,
grid = :on,
left_margin=2cm, right_margin=8cm, top_margin=2cm)
p = twinx()
plot!(p, a, log.(a),
label = "log(x)",
ylabel = "The right Y label",
color = :green,
legend = :topright,
tick_direction = :out,
# grid = :on,
left_margin=2cm, right_margin=8cm, top_margin=2cm,
box = :on)
savefig("tstplt6.png")
endfunction tstplt7()
a = 1:3
b = [1, 2, 7]
plot(a, b,
label = "randData",
xlabel = "numbers",
ylabel = "Rand data",
color = :red,
legend = :topleft,
grid = :on, left_margin=2cm, right_margin=8cm, top_margin=2cm)
plot!(twinx(), a, log.(a),
label = "log(x)",
ylabel = "The right Y label",
color = :green,
legend = :bottomright,
grid = :on, left_margin=2cm, right_margin=8cm, top_margin=2cm)
savefig("tstplt7.png")
end# bok: Body Of Knowledge
# iPoint: index specifying last point to plot
# example: bok()
function bok(iPoint::Int64=13)
# initialize
n = 12
xmax = div(n^2,10,RoundUp) * 10
ymax = div(4*n+2,10,RoundUp) * 10
plot(seriestype = :scatter,
xaxis = :log10,
yaxis = :log10,
xlims = (0.4, 180),
ylims = (0.8, 120),
xticks = ([0.5, 1, 2, 10, 100, n^2],
[latexstring(L"-\infty"),
latexstring(L"10^0"),
"2",
latexstring(L"10^1"),
latexstring(L"10^2"),
latexstring(L"n^2")]),
xlabel = latexstring("common ratio ", L"r", "\ndegrees of freedom"),
ylabel = latexstring("coefficient ", L"a", "\ndegrees of freedom"),
grid = :on,
left_margin=2cm, right_margin=8cm, top_margin=2cm)
plot!(twinx(),
seriestype = :scatter,
xaxis = :log10,
yaxis = :log10,
xlims = (0.4, 180),
ylims = (0.8, 120),
yticks = ([1, 6, 8, 10, n+1, 4*n+2, 100],
["",
"6",
"8",
"",
latexstring(L"n+1"),
latexstring(L"4n+2"),
""]),
grid = :on,
left_margin=2cm, right_margin=8cm, top_margin=2cm)# geometric series
A = [1]
R = [1]
plot!(R, A,
marker = :circle,
markersize = 16,
color = :white,
markerstrokecolor = :red,
label = "geometric series")
# geometric series closed form
A = [1]
R = [1]
plot!(R, A,
xaxis = :log10,
yaxis = :log10,
marker = :utriangle,
markersize = 10,
color = :red,
linecolor = :white,
label = "geometric series closed form")
if (iPoint <= 2)
@goto printplot
end
# complex geometric series
A = [1]
R = [2]
plot!(R, A,
marker = :circle,
markersize = 16,
color = :white,
markerstrokecolor = :green,
label = "complex geometric series")
# complex geometric series closed form
A = [1]
R = [2]
plot!(R, A,
marker = :utriangle,
markersize = 10,
color = :green,
linecolor = :white,
label = "complex geometric series closed form")
if (iPoint <= 4)
@goto printplot
end
# Laurent series
A = [4*n+2]
R = [2]
plot!(R, A,
marker = :circle,
markersize = 16,
color = :white,
markerstrokecolor = :red,
label = "Laurent series")
# complex Fourier series
A = [4*n+2]
R = [2]
plot!(R, A,
marker = :pentagon,
markersize = 10,
color = :red,
linecolor = :white,
label = "complex Fourier series")
if (iPoint <= 6)
@goto printplot
end
# power series
A = [n+1]
R = [1]
plot!(R, A,
marker = :circle,
markersize = 16,
color = :white,
markerstrokecolor = :blue,
label = "power series")
# Taylor series
A = [n+1]
R = [1]
plot!(R, A,
marker = :square,
markersize = 10,
color = :blue,
linecolor = :white,
label = "Taylor series")
if (iPoint <= 8)
@goto printplot
end
# binary encoded numbers
A = [n+1]
R = [0.5]
plot!(R, A,
marker = :circle,
markersize = 12,
color = :white,
markerstrokecolor = :black,
markerstrokewidth = 3,
label = "binary encoded numbers")
annotate!(0.65, 0.8, # impose text white background
text("█",
:white,
:center))
annotate!(0.65, 0.8, # denote axis split
text("/ /",
:black,
:center))
if (iPoint <= 9)
@goto printplot
end# Bezier curves
A = [6]
R = [1]
plot!(R, A,
marker = :utriangle,
markersize = 12,
color = :white,
markerstrokecolor = :black,
markerstrokewidth = 3,
label = "quadratic Bezier curves")
A = [8]
R = [1]
plot!(R, A,
marker = :square,
markersize = 12,
color = :white,
markerstrokecolor = :black,
markerstrokewidth = 3,
label = "cubic Bezier curves")
if (iPoint <= 11)
@goto printplot
end
# matrix polynomials
A = [1]
R = [n^2]
plot!(R, A,
marker = :circle,
markersize = 16,
color = :white,
markerstrokecolor = :green,
label = "matrix polynomials")
# matrix exponentials
A = [1]
R = [n^2]
plot!(R, A,
marker = :hexagon,
markersize = 10,
color = :green,
linecolor = :white,
label = "matrix exponentials")
@label printplot
if (iPoint < 13)
fname = @sprintf("bok%02d.png", iPoint)
else
fname = "bok.png"
end
savefig(fname)
endrectangle(w, h, x, y) = Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])function fp32()
# initialize
nBit = 32
x0 = 1
x1 = 15
y0 = 3
y1 = 4
plot(legend = false,
showaxis = false,
size = (800, 300),
ticks = false,
xlims = (0, 16),
ylims = (0, 6),
title = "binary32 IEEE 754 standard\n" *
"single-precision floating point number format")
# color bit fields
plot!(rectangle(23*(x1-x0)/nBit,1, x1-23*(x1-x0)/nBit,y0),
opacity = 0.3,
color = :red)
plot!(rectangle(8*(x1-x0)/nBit,1, x0+(x1-x0)/nBit,y0),
opacity = 0.3,
color = :green)
plot!(rectangle((x1-x0)/nBit,1, x0,y0),
opacity = 0.3,
color = :blue)
# draw bit boundaries
plot!([x0; x1], [y0; y0],
color = :black)
plot!([x0; x1], [y1; y1],
color = :black)
for iBit=0:nBit
x = x1 - iBit*(x1-x0)/nBit
plot!([x; x], [y0; y1], color = :black)
if (iBit < nBit)
annotate!(x-(x1-x0)/(2*nBit), y1+0.5,
text(string(mod(iBit,10)), :black, :hcenter))
end
if (mod(iBit,10) == 0)
annotate!(x-(x1-x0)/(2*nBit), y1+1.2,
text(string(div(iBit,10)), :black, :hcenter))
end
end
# write labels
annotate!(x0, y1+0.8, text("bit #:", :black, :right))
annotate!(x1-23/2*(x1-x0)/nBit, y0-0.5, text("fraction (23 bits)", :black, :hcenter))
annotate!(x0+(1+8/2)*(x1-x0)/nBit, y0-0.5, text("exponent (8 bits)", :black, :hcenter))
annotate!(x0+(x1-x0)/(2*nBit), y0-0.7, text("sign", :center, rotation = 90))
fname = "fp32.png"
savefig(fname)
endtrapezoid(x0, y0, x1, y1) = Shape([x0,x0,x1,x1], [0,y0,y1,0])function selfsim()
print("selfsim")
# initialize
CSEQ = [:red, :green, :blue] # Color SEQuence
nTrap::Int32 = 25
nFrame::Int32 = 80
xmaxmax = 2.0
xmaxmin = 1.0
dgrid = 0.5
r = (1/2)^(2/3)
# for each animation frame
for iFrame = 0:nFrame-1
# display progress
if (mod(iFrame+1,10) == 0)
print(".")
end
# start plot of new frame
xmax = xmaxmax - (xmaxmax-xmaxmin)*iFrame/nFrame
# print("\n", xmax, ">")
plot(legend = false,
box = true,
ticks = false,
aspect_ratio = :equal,
xlims = (0, xmax),
ylims = (0, 1.5*xmax))
# prepare first (smallest) trapezoid
x0 = 2*r^(nTrap/2)
y0 = r^(nTrap/2)/(1-r)
x1 = 2*r^((nTrap-1)/2)
y1 = r^((nTrap-1)/2)/(1-r)
# for each trapezoid
for iTrap = nTrap-1:-1:0
iColor::Int16 = mod(iTrap,3) + 1;
# color trapezoid area
plot!(trapezoid(x0,y0, x1,y1),
opacity = 0.3,
color = CSEQ[iColor])
# prepare next trapezoid
x0 = x1
y0 = y1
x1 = 2*r^((iTrap-1)/2)
y1 = r^((iTrap-1)/2)/(1-r)
end
fname = @sprintf("./selfsim/frame%05d.png", iFrame);
savefig(fname)
end
print("\n")
# /tool/ffmpeg-4.3.1/bin/ffmpeg -i _selfsim.mp4 -r 12 _selfsim.gif
# /tool/ffmpeg-4.3.1/bin/ffmpeg -f image2 -i frame%05d.png -vf scale=1080:-2,setsar=1:1 -f mp4 -q:v 0 -vcodec mpeg4 -r 20 _selfsim.mp4
# see https://trac.ffmpeg.org/wiki/Encode/VP9
end# tsa: Taylor Series Approximation
# nDegree: degree of last term in partial series
# example: tsa()
function tsa(nDegree::Int64=19)
# define coefficients (1-based index)
C = [
0, 1/factorial(1),
0, -1/factorial(3),
0, 1/factorial(5),
0, -1/factorial(7),
0, 1/factorial(9),
0, -1/factorial(11),
0, 1/factorial(13),
0, -1/factorial(15),
0, 1/factorial(17),
0, -1/factorial(19)]
if (nDegree >= length(C))
@printf("error: nDegree cannot be more than %d\n",
length(C)-1)
return
end
# define sin function
nPoint = 201
X = LinRange(-1.0, 1.0, nPoint)
F = sin.(X)
# define Taylor series approximation
TSF = zeros(nPoint, nDegree+1)
for iDegree = 0:nDegree
TSF[:,iDegree+1] = C[iDegree+1] .* X .^ iDegree
end
TSFC = cumsum(TSF, dims=2)
# show Taylor series approximation error
p1 = plot(X, F,
legend = false)
p2 = plot(X, TSFC[:,end] - F,
legend = false)
M = maximum(abs.(TSFC-repeat(F,1,nDegree+1)), dims=1)
display(M)
p3 = plot(M',
legend = false,
yaxis = :log)
plot(p1, p2, p3, layout=(3,1))
end
