Math Programming

Creation of MathGraph3D (Part 1 – Foundation)

This series focuses on the creation of the original version of my 3D plotting software MathGraph3D. This first part is concerned with the overall structural components of the software.

Before getting into any actual algorithms or math, it’s necessary to set up an outline for how the code will be structured. I decided to create a manager class called “Plot” that would handle everything to do with 3D space and all objects in the calculator. There are five main operations that the Plot must be able to do (among other smaller operations):

  • Determine the bounds and scale of the 3D space. The user should be able to input upper and lower bounds for each axis, and this defines the space used by the calculator.
  • Projects points in 3D to cartesian points in 2D, then finally to canvas coordinates, relative to the top left of the screen, to be drawn. Typically, this is done with a camera and matrices for rotation and projection. But the original version of MathGraph did not use this because I wasn’t familiar with matrix math, so I developed my own algorithm for projection. However, both will be covered in this series for the sake of its usefulness.
  • Occlusion of shapes. I decided to go with the Painter’s method for occlusion. Since MathGraph3D does not operate with pixel-perfect access, better algorithms are mostly out of reach. “Occlusion” means that when an object (polygons and lines) is obscured or blocked by another object from the perspective of the view, it is not drawn to the screen. However, the Painter’s method is imperfect and sometimes results in ambiguous situations. Take, for example, the scenario in which two polygons intersect at a line through their centers. There is no way to decide which polygon should be in front because both are partially obscured by the other. To handle this, I developed a modification on the standard algorithm that will be discussed later.
  • Drawing Everything to the Canvas. This includes functions, points, vectors, geometry, the axes, axis numbers, etc.
  • Clipping. Objects on the Plot should never exceed the bounds of the defined 3D space. For example, assume all axes range from -4 to 4. The function $f(x,y)=x^2+y^2$ will exceed the upper z bound at input values outside the circle centered at the origin with radius 2. But I don’t want the part of the surface that goes outside the bounds to be shown; it should be “clipped” at the boundary. The same goes for all bounding planes (x=-4, x=4, y=-4, y=4, z=-4, z=4) and any additional ones defined by the user.

The next step for baseline setup is the general “Plot object” class, from which every surface should inherit. This will be covered in the next post. The final main structural component of the code is the GUI, but that is not discussed here. Now, we can get into the details of the Plot object. I’ll discuss the main operations of the PLot in the order they are outlined above.

Determining the bounds and scale of the 3D space

This operation involves determining the position of the axes and the origin based on bounds, the scale per unit in each direction based on the axis length, and the number of units in each direction.

It is easy to find the number of units; it’s just the length of the interval between the upper and lower boundaries. To compute that, you can subtract the lower boundary from the upper one. To find the unit scales, you can take the axis length and divide it by the number of units along that axis. For example, if one of the axis lengths is 500 pixels and the user wants it to range from -5 to 5, there are 5 – (-5) = 10 units, so there are 500 pixels / 10 units = 50 pixels per unit.

The problem is actually a bit more complicated than that because all axes are not necessarily the same length but they should be the correct length relative to each other. In other words, if one axis has 10 units and another has 5, the first should appear twice as large as the first. Solving this will automatically yield the correct position of the origin.

3D Projection – version implemented in the original MathGraph3D

I’ve written about my algorithm before, and it’s buried deep somewhere on this website as a PDF. It’s not very high quality, so I’ll take this opportunity to rewrite it more clearly.

Figure 1 shows the three dimensional axes $x_3$, $y_3$, and $z_3$. These are the axes for the three dimensional points to be projected. $\hat{i}, \hat{j}, \hat{k}$ are the units of the $x_3$, $y_3$, and $z_3$ directions, respectively. $x_2$ and $y_2$ are the two-dimensional axes onto which the 3D points are to be projected.

Since any point can be expressed as a linear combination (scaled sum) of the unit vectors, the problem of projecting any 3D point to 2D is a matter of finding the projected unit vectors.Then any 3D point $\left(a, b, c\right)$ may be projected like this, where $\text{P}=\left(x_p, y_p\right)$ is the result: 

$$x_p=a\cdot \hat{i}_x + b \cdot \hat{j}_x + c\cdot \hat{k}_x \\y_p=a\cdot\hat{i}_y + b\cdot\hat{j}_y + c\cdot\hat{k}_y$$

Written in matrix form,

$$\text{P}=\begin{bmatrix}a&b&c\end{bmatrix}\begin{bmatrix} \hat{i}_x & \hat{i}_y \\ \hat{j}_x & \hat{j}_y \\ \hat{k}_x & \hat{k}_y \end{bmatrix}=\begin{bmatrix}a\hat{i}_x+b\hat{j}_x, a\hat{i}_y+b\hat{j}_y+c\hat{k}_y\end{bmatrix}$$

It’s worth noting that in the original version’s code, the multiplication of $c$ and $\hat{k}_x$ is not actually computed because the z axis is fixed to always be straight up and down. This means $\hat{k}_x$ is always going to be zero so doing this computation is a waste.

Let $y’$ be the line normal to the plane formed by $x_2$ and $y_2$ and let $\alpha$ be the angle between $x_2$ and $y_2$ in the plane formed by $x_2$ and $y’$. Then $\alpha$ corresponds to rotation about the z axis in 3D space. Let $\beta$ be the angle formed between $y_2$ and $z_3$ in the plane formed by $y_2$ and $y’$. So $\beta$ corresponds to rotation about a stationary horizontal axis in the 3D space.

$\hat{k}$ is the easiest unit vector to find because its projected x value never changes.

In Figure 2 it can be seen that:

$$\hat{k}_y = \cos\left(\beta \right) \implies \hat{k}_{\text{projected}} = \left(0, \cos\left(\beta \right) \right)$$

Finding the projected x of $\hat{i}$ is straightforward, but the projected y is a bit more difficult.

Figure 3 shows that $\hat{i}_x=\cos\left(\alpha \right)$. If we take the side of the triangle with $\sin\left(\alpha \right)$ and project it onto the plane formed by $y_2$ and $y’$ and angle $\beta$, this triangle is obtained:

Figure 4 shows:

$$\sin\left(\beta \right) = \frac{\hat{i}_y}{\sin\left(\alpha \right)} \implies \hat{i}_y = \sin\left(\alpha \right)\sin\left(\beta \right)$$

$$\text{So } \hat{i}_{\text{projected}} = \left(\cos\left(\alpha \right),\sin\left(\alpha \right)\sin\left(\beta \right)\right)$$

The angle in the same plane as $\alpha$ but between $x_2$ and $y_3$ is $\alpha+\frac{\pi}{2}$ because $x_3$ and $y_3$ are perpendicular. Therefore $\hat{j}_x$ must be

$$\cos\left(\alpha + \frac{\pi}{2}\right) = -\sin\left(\alpha \right)$$

Since $\sin\left(x+\frac{\pi}{2}\right)=\cos\left(x\right)$. So the same triangle as in figure 4 can be used except with $\cos\left(\alpha \right)$ instead of $\sin\left(\alpha \right)$ and $\hat{j}_y$ instead of $\hat{i}_y$:

We can see from Figure 5 that:

$$\sin\left(\beta \right) = \frac{\hat{j}_y}{\cos\left(\alpha \right)} \implies \hat{j}_y = \cos\left(\alpha \right)\sin\left(\beta \right)$$

$$\text{So } \hat{j}_{\text{projected}} = \left(-\sin\left(\alpha \right), \cos\left(\alpha \right)\sin\left(\beta \right) \right)$$


$$\hat{i}_{\text{projected}} = \left(\cos\left(\alpha \right),\sin\left(\alpha \right)\sin\left(\beta \right)\right)$$

$$\hat{j}_{\text{projected}} = \left(-\sin\left(\alpha \right), \cos\left(\alpha \right)\sin\left(\beta \right) \right)$$

$$\hat{k}_{\text{projected}} = \left(0, \cos\left(\beta \right) \right)$$

Occlusion of shapes

In order to implement the Painter’s method, we must have a way to determine which shapes are closest to the view at any angle. Initially, I went case by case and determined the closest octant to the viewer and sorted the polygons according to its outer corner. This, however, is imperfect. So once again I began to create my own algorithm. I’m somewhat of a chaotic worker with most of my progress and ideas coming in bursts, which means when inspiration strikes I have to find the quickest way to get the idea out on paper before I forget it. Unfortunately, I can’t find the scrap of paper where I originally derived my algorithm for this, and at the moment I can’t rederive it, so I’ll present the formula here and hopefully update this post if I find the paper or remember my train of thought.

The idea behind this algorithm is to find the point of intersection between a ray originating at the viewer’s eye and a sphere centered at the origin with given radius $r$. The polygons can then be sorted in reverse by squared distance from the intersection point and they will be drawn in the correct order, accounting for occlusion. For efficiency reasons, squared distance is used instead of actual distance, as the square root calculation is computationally expensive. The squared distance will preserve the correct order, as $\sqrt{x}$ is injective and monotonically increasing across the positive numbers. Anyway, suppose the intersection point is called $Q$. Let $\theta = -\left(\alpha +\frac{\pi}{2}\right)$ and $r$ be the radius of the sphere centered at the origin. Then the coordinates of $Q$ are given by:

$$Q = \left( r\cos\left(\theta \right), r\sin\left(\theta \right), \left| r \right| \sin\left(\pi – \beta \right) \right)$$

I’m not sure what the absolute value signs are there for, because this formula works fine for any positive radius. It’s probably there in case I wanted to have a negative radius for some reason, but it’s perfectly acceptable to leave it out.

The modification on the Painter’s algorithm that I referred to earlier is splitting polygons to handle overlap. This modification was not included in the first version, but it’s worth including here. It works by checking every single polygon against every other one, and if there’s overlap, splitting both polygons into two pieces at the line of intersection. Solving for the line of intersection between two planes in space is a basic linear algebra problem, and there are many resources online that show how to compute it so I won’t cover it here. The only difference is that the polygons are plane regions, not entire planes. But resolving this difference is only a matter of restricting the boundaries of the plane on which the polygon lies.

Surface Clipping

The original version of MathGraph3D only allowed clipping for planes equal to z=c for some constant c. I’ll cover the general plane clipping method in this post anyway. We want an algorithm that can take in the normal vector of a plane and a point that lies on the plane, and replace all points on the side to which the normal points (the “outside”) with interpolated points on the plane itself. This is done by first identifying which points are on the outside, then computing the intersection of the segment formed between each outside point and the nearest inside point of a polygon. If all points are outside, the polygon is removed entirely; if no points are outside, the polygon is kept as is. Call the normal vector $\vec{n}$ and the point lying on the plane $\vec{R}$. Then, the clipping plane precomputes the test value $T = \vec{n} \cdot \vec{R}$. All points on the outside of the plane, when dotted with the normal vector, will produce a value greater than $T$. Therefore, the following test can be used for a point $\vec{S}$:

$$\vec{n}\cdot \vec{S} – T > 0$$

If the test is true, $S$ is outside the clipping plane. Otherwise, it is on the inside. Now, given an outside point $P_{\text{out}}$ and an inside point $P_{\text{in}}$, we need to find the point on the clipping plane that lies between them. First, the direction vector of the line between these points is computed: $\vec{d} = P_{\text{out}} – P_{\text{in}}$. We then compute the value of the parameter $t$ at which the line between these two points intersects the clipping plane: 

$$t=\dfrac{\vec{n} \cdot (\vec{R} – P_{\text{out}})}{\vec{n}\cdot \vec{d}}$$

From here, computation of the intersection point is easy:

$$ P_{\text{out}} + t\cdot \vec{d} $$

Each new polygon is made up of the inside points and the computed boundary points, all sorted clockwise.


Building a 3D Boat in Desmos

The graphing calculator Desmos recently held a competition called the Global Math Art Contest. People aged 13-18 submitted artwork they created with the graphing calculator using combinations of math curves. The judging is based on creativity, originality, visual appeal, and the math used to create the art. As a person with little artistic talent, I had to rely on the use of math category to get my edge. So I decided to create a rotatable 3D boat that moves and tilts according to the waves in the water below it. Here I’ll explain how I went about doing this. Finished project:

“Modifying” the in-browser calculator

By default, Desmos only supports 6 different colors and none of them are really optimal for a boat. I wanted to use brown for the boat and light blue for the water beneath it. I opened developer tools in Chrome to see how I could add additional colors. As it turns out, Desmos defines an instance of the calculator called Calc that can be changed from the client side. After looking at the Desmos API documentation for a little bit I was able to figure out how to modify the available colors. I used a Chrome extension called Tampermonkey to automatically execute my color adding script after the calculator loads. Here’s the script:

Doing this might actually disqualify me from the contest.. but I checked the rules and there was no mention of using the Desmos API. It just said that the art had to be done in the online calculator, which it was. I won’t know for sure if it’s allowed until the results are announced in the second half of May.

Preparing to draw in 3D

Since Desmos is strictly a 2D graphing calculator, I had to implement my own algorithms for drawing in three dimensions. I’ve already worked out how to project from 3D world space to a 2D canvas on a screen because I needed to do that for MathGraph3D. The algorithm I implemented in Desmos is actually the OLD version of my 3D to 2D algorithm. Here is the derivation.

The basic idea is to figure out the screen coordinates of each 3D unit vector. First, we must have a way to define orientation. Call $\alpha$ the angle formed between the positive $x$ axis and a fixed vertical axis, and call $\beta$ the angle between the positive $x$ axis and a fixed horizontal axis. If $\vec{i}$, $\vec{j}$, and $\vec{k}$ are the standard unit vectors for $\mathbb{R}^3$, their projections are

$$ \vec{i}_{\text{proj}} = \left(\cos\left(\alpha\right),\sin\left(\alpha\right)\sin\left(\beta\right)\right)$$

$$ \vec{j}_{\text{proj}} = \left(-\sin\left(\alpha\right),\cos\left(\alpha\right)\sin\left(\beta\right)\right) $$

$$\vec{k}_{\text{proj}} = \left(0,\cos\left(\beta\right)\right) $$

Now, the boundaries of the 3D space must be defined. To keep things simple I decided to make the space a cube. So in the graph the value $s_{tart}$ is necessarily negative and defines the lower bound for the space in each dimension. Then $s_{top} = -s_{tart}$ defines the upper bound and $u_{nits} = s_{top} – s_{tart}$ is the number of units in each direction. I disabled the gridlines, axis numbers, and axes so the graph is just a blank white page by this point.

The water

The top of the water is the surface of a function of two variables. I wanted it to look quite realistic, so I experimented with some functions to create ripples and waves. I took advantage of the wavy nature of sinusoid functions to do this. But that comes with a different challenge: sines and cosines are periodic, and I didn’t want the waves to look like a simple pattern that repeats itself over and over. To avoid this, I messed with the period and amplitude of a sine wave and a cosine wave, then added a constant to the argument of the cosine to translate it in and out of phase with the sine. When the constant’s value is animated, the function looks like this:

(This function is defined by $w_{1}\left(x\right)=\frac{1}{4}\left(\sin\left(x\right)+\frac{1}{2}\cos\left(2x-b_{0}\right)\right)$ where $b_0$ is the constant)

Next, I defined another function $w_{2}\left(x\right)$ that is exactly the same as $w_{1}\left(x\right)$ except the constant $b_0$ is subtracted from the argument of cosine, not added. Finally, to make an occasional really big wave, I defined another function $w_{3}\left(x\right)=2e^{-\frac{1}{16}\left(x-a_{0}\right)^{2}}$. This is a scaled normal distribution with mean $a_0$. When $a_0$ is animated, the function looks like this:

These are all the tools needed to create the final water function $W\left(x,y\right)$. Here is its definition:

$$|x+y|+|x-y|\le u_{nits}:w_{1}\left(x\right)+w_{2}\left(y\right)+w_{3}\left(x+y\right)$$

The condition on that function tests if the point $(x,y)$ is in the portion of the $xy$-plane being drawn to the screen. This is needed for drawing the polygons at the side of the water to shade from the surface of the water to the bottom of the space. Notice that $w_1$ acts in the $x$ direction, $w_2$ acts in the $y$ direction, and $w_3$ acts in the $x+y$ direction, or diagonally. This creates more wave interference and thus further removes the feeling that the ripples are just a repeated pattern. Here’s what it looks like so far:

Change of basis

In reality, when I was making this project I created the boat first then added change of basis later. However, it makes more sense to explain it this way. Change of basis in this context means defining a new linearly independent set of basis unit vectors for $\mathbb{R}^3$ such that the $z$ unit vector is normal to the surface of the water under the boat. Let the new basis be $\{e_1,e_2,e_3\}$. Then $e_1$ is calculated by taking the vector difference between $(-l, 0, W(-l,0))$ and $(l, 0, W(l,0))$ where $l$ is the length of the boat. $e_2$ is the vector difference between $(0, -w, W(0, -w))$ and $(0, w, W(0, w))$ where $w$ is the width of the boat. $e_3$ is the cross product $e_1 \times e_2$. Finally, these vectors are all normalized.

Now, by multiplying the coordinates of each point on the boat by the corresponding vector in the new basis, the boat will tilt in both directions according to the shape of the water beneath it. It’s a bit hard to tell that it’s working until the big wave comes by and the boat dramatically tilts forward.

The boat

This was by far the hardest part because I’m not an artist and ruining the project with a poorly done boat would have been disappointing. But I managed to put together a pretty decent boat that, honestly, turned out better than I expected. I defined three functions $g_1$, $g_2$, and $g_3$. $g_1$ defines the shape of the boat along the front to back direction, $g_2$ defines the shape along the side to side direction, and $g_3$ is used to make the boat have curved (or “flared-out”) edges instead of straight sides. $g_1$ is a parabolic segment that is divided by the square of the boat length, so that it’s slightly curved upward but not too much. It looks like this:

$g_2$ is a hyperbolic cosine function divided by the square of the boat width. I tried a parabola at first but it was too round near the center.

$g_3$ is another normal distribution curve that depends on the width of the boat:

Putting these together, and using Desmos list magic to draw many of these curves, the boat looks like this:

Finally, I used a grid of sample $W$ values to determine the maximum height of the water near the boat. I translated the boat upwards by this amount divided by two so it stays near the surface but looks like it is displacing water.

The link to the final graph is near the top of this post. Thanks for reading!


Common misconceptions about the differentiation operator

(Note: math may not render properly in Microsoft Edge. Any other browser should work, though. I’m trying to figure out the source of this issue.) This is the differentiation operator in single variable calculus:


Since several procedures in calculus happen to involve manipulating the “numerator” and “denominator” of the differentiation operator, people have gained the incorrect idea that it IS in fact a ratio of two “infinitely small” quantities. This is absolutely not true! Modern mathematics does indeed have a notion of infinitely small quantities (called surreal numbers), but they are defined within an entirely different mathematical system called nonstandard analysis. The whole motive behind formulating a rigorous definition of a limit is avoiding problems with infinity. That is to say, $\dfrac{dy}{dx}$ is not a ratio of an infinitely small change in $y$ to an infinitely small change in $x$; instead, it is shorthand for the limit of the difference quotient of some function $y(x)$. By definition,

$$\frac{d}{dx}{y(x)}=\lim_{\Delta x\to 0}\frac{y(x+\Delta x) – y(x)}{\Delta x}$$

The limit denotes the quantity that the difference quotient “approaches” as $\Delta x$ grows small; that is, for any arbitrarily small difference $\varepsilon$ between the value of $\frac{y(x+\Delta x)-y(x)}{\Delta x}$ and the limit $\lim_{\Delta x\to 0}\frac{y(x+\Delta x) – y(x)}{\Delta x}$, there exists a choice of $\Delta x$ that yields a value of the difference quotient within $\varepsilon$ of the limit. Nothing in this definition has anything to do with “infinitely small quantities”. Instead, it deals with “arbitrarily small” quantities. In the end, though, the limit is a different number altogether that just so happens to be associated with the expression. The expression can get within arbitrarily small distance of the limit given an appropriate choice of $\Delta x$.

Now that the definition is out of the way, it seems necessary to address several instances in calculus where it seems okay to manipulate the differentiation operator as if it’s a fraction. I will cover the the most common situation (and ALL others can be explained, rigorously, using a concept called “differential forms”, which is currently beyond my understanding). The two most common are in the separation of variables for solving simple differential equations and the method of integration by substitution (which I was planning to cover as well, but it turns out Wikipedia has an excellent explanation here).

Separation of variables

After separating the variables in a differential equation, we are left with something of the form:


At which point we can take the antiderivative of both sides with respect to $x$ (since antiderivatives are unique up to a constant):


The right hand side of this equation essentially asks, “what function, when we take its derivative with respect to $x$, yields $f(y)\frac{dy}{dx}$”. If we let


then, by the chain rule,


which is of course the equivalent to


Therefore, integrating both sides with respect to $x$,


which justifies the final equation


So in practice, it seems like you’re “multiplying both sides of the equation by $dx$”, but that is the consequence of the work shown above.

Math Programming

Fundamental Theorem of Calculus Explained

Relation between derivatives and integrals

Why are derivatives the inverse of integrals?

By Sam Brunacini

Understandably, most people get confused when first introduced to the idea that finding the area under a curve is the inverse to finding instantaneous rate of change. It seems like the intuition behind this is hidden behind all the formulas they make you learn. Here I’ll explain why the Fundamental Theorem of Calculus (FTC) works. This is not a rigorous proof of the theorem, just an explanation of why they make sense.

First, let’s write the FTC: $$\text{Part 1: } \frac{d}{dx}\int_{a}^{x}f(t)dt = f(x)\\$$ $$\text{Part 2: } \int_{a}^{b}f(x)dx=F(b)-F(a)$$

To help understand Part 1, imagine you own a library. Before the library opens at 8:00 AM, you count 100 books in stock. Between 8:00 and 9:00, 5 books are withdrawn. Between 9:00 and 10:00, 8 books are returned. Without recounting, what is the number of books in stock at 10:00? Of course, you subtract 5 from 100 to account for the withdrawals then add 8 for the returns. This leaves $100 – 5 + 8 = 103$.

Believe it or not, this is exactly what FTC Part 1 says. Call the number of books in the library $x$ hours after 8:00 $f(x)$. Then $f(0)=100$ because that’s how many books there were before anyone came in. $f(x)$ changes by -5 during the first hour. We then add this change to the total to get $f(1)=95$. Due to the returns, $f(x)$ changes by 8 in the second hour. Adding this change to the total, we know $f(2)=103$. FTC part 1 says that a function is equal to the accumulation (or total) of the changes in itself as $x$ varies. The derivative indicates change, and the integral indicates accumulation, or summing the values together. This is exactly what happened in the library example.

In [47]:
import matplotlib.pyplot as plt
import numpy as np
import math

THIRD_PI = math.pi / 3

plt.rcParams['figure.figsize'] = [15, 5]
fig, (ax1, ax2, ax3) = plt.subplots(1, 3)
fig.suptitle("The finite accumulation (sum) of changes")
ax1.set_title("Sum of 4 changes")
ax2.set_title("Sum of 6 changes")
ax3.set_title("Sum of 10 changes")

def plot_cos(axes, steps=100, color="blue"):
    X = np.linspace(0, THIRD_PI, num=steps)
    Y = np.cos(5*X)
    axes.plot(X, Y, color=color)

plot_cos(ax1, color="red")
plot_cos(ax2, color="red")
plot_cos(ax3, color="red")

plot_cos(ax1, 4)
plot_cos(ax2, 6)
plot_cos(ax3, 10)

The graphs above show finite sums of changes in some function. Roughly speaking, an integral is the sum of infinitely many tiny changes in the function. Notice that as the number of changes increases, the approximation (blue line) gets closer to the exact curve (red line). In fact, the approximation gets arbitrarily close to the curve as the number of changes grows larger.

Now, for Part 2. Remember that $F(x)$ is defined as the function that gives the area under a curve from 0 to $x$. Assume for now that $a \lt b$. This is a safe assumption since we can use the formula $$\int_{b}^{a}f(x)dx=-\int_{a}^{b}f(x)dx$$ to make it true if originally $a \gt b$. $F(a)$ is then the area under the curve from 0 to $a$ and $F(b)$ is the area under the curve from 0 to $b$. If we plot these separately, it looks like this:

In [55]:
import matplotlib.pyplot as plt
import numpy as np
import math

THIRD_PI = math.pi / 3

plt.rcParams['figure.figsize'] = [15, 5]
fig, (ax1, ax2) = plt.subplots(1, 2)

def plot_sin(axes, steps=100, color="red"):
    X = np.linspace(0, THIRD_PI, num=steps)
    Y = np.sin(4*X)/2+0.5
    axes.plot(X, Y, color=color)

X1 = np.linspace(0, THIRD_PI/3)
ax1.fill_between(X1, np.sin(4*X1)/2+0.5, color="blue", alpha=0.25)
ax1.annotate("x=a", xy=(THIRD_PI / 3, math.sin(4 * THIRD_PI / 3) / 2 + 0.4), size=15)

X2 = np.linspace(0, THIRD_PI*0.75)
ax2.fill_between(X2, np.sin(4*X2)/2+0.5, color="green", alpha=0.25)
ax2.annotate("x=b", xy=(THIRD_PI * 0.75, math.sin(3*THIRD_PI)/2+0.5), size=15)

If you subtract the area shown in the first diagram from the second (which is simply doing $F(b) – F(a)$), you end up with the area between $a$ and $b$.

In [ ]:
Math Programming

Fluid flow in vector fields

I made a python program that simulates fluid particles flowing through a vector field. Here I’ll put some things I discovered, some example gifs, and the complete code file.

Some patterns

A vector field is created by inputting points into a vector valued function and placing the tail of each output vector at its corresponding input. As I explored the particle flow of vector fields of different functions with my program, a few patterns started to emerge. The main ones are:

  • Attraction points – points on the vector field where particles tend to converge (I say “tend to” because the particles don’t always flow directly towards the attraction points, but as time goes on, they get closer).
  • Attraction lines – entire curves made up of attraction points
  • Repelling points – points where particles tend to diverge
  • Repelling lines – entire curves made up of repelling points
  • Swirly points – regions where the fluid particles swirl around but do not get pulled in.
  • I bet there’s such a thing as swirly lines also. But I haven’t been able to think of a function where this occurs.

(Note that I made these names up – I haven’t ever studied vector fields formally)

Here’s an example of each of these patterns.

Attraction points – $\vec{f}(x, y)=\left \langle \arctan(y-x), x\cos(x) \right \rangle$

Attraction curve – $\vec{f}(x, y)=\left \langle \cos(x^2+y^2),\sin(x^2+y^2) \right \rangle$

Fluid flow 2

Repelling point – $\vec{f}(x, y)=\left \langle x-y, x+y \right \rangle$

Fluid flow 3

Repelling Curve – $\vec{f}(x, y)=\left \langle 0, y+\sin(x) \right \rangle$

Fluid flow 4

Both repelling points and attraction points – $\vec{f}(x, y)=\left \langle \sin(x-y), \sin(x+y) \right \rangle$

Fluid flow 5

Swirly point – $\vec{f}(x, y)=\left \langle \sin(2y), -\frac{x}{2} \right \rangle$

Fluid flow 6

The code

Here’s the complete code:


Collection of Desmos Graphs

Recently I went through my old Desmos ( graphs and I’d like to share 5 of my favorites here. There are some additional interesting graphs at the bottom that are cool as well. You can click on the links after each header to edit and interact with the graphs in the browser.

Complex Plane Rotation

Try it yourself:

The goal of this graph is to take a function $f(x)$ and rotate it around the origin by an arbitrary radian amount in order to demonstrate the concept of rotation in the complex plane. I accomplished this by first looking to rotations in the complex plane, then translating the method to the real plane.

The blue curve is the original function and the orange curves are the rotated versions.

A complex number $z=x+yi$ is plotted in the complex plane by putting the real part $x$ on the horizontal axis and the imaginary part $y$ on the vertical axis. Note that this means each “point” is really just a single number. This means a $z$ can be rotated about the origin $a$ radians by multiplying by $i^{\frac{2a}{\pi}}$. For example, consider $z=1+0i=1$. If you let $a=\pi /2$, the factor becomes $i$, which corresponds to a ninety degree rotation. If you let $a=\pi$, the factor becomes, $i^{2}=-1$, which corresponds to a 180 degree rotation. Complex numbers can also be represented in polar form by $z=r\cos{\left(\theta\right)}+r\sin{\left(\theta\right)}i$, where $r=\sqrt{x^2+y^2}$ (the magnitude of the vector $x\hat{\imath}+y\hat{\jmath}$) and $\theta$ is the angle of the same vector with the positive x axis.

Unfortunately, Desmos does not support complex numbers. Luckily, the parameterization of the complex unit circle is exactly the same as the parameterization of the unit circle in the real plane (excluding the imaginary unit of course). In the Desmos graph I took advantage of this fact to rotate each point. To rotate an entire curve (as opposed to a single point), I made a parametric function of $t$ and applied the transformation to each point $(t, f(t))$ in the domain $a\leq t \leq b$.

Tangent and normal lines

Try it yourself:

I like this one because of the simplicity. I also like it because I had no idea how to do calculus when I created it about a year and a half ago. However, I did understand that the derivative magically yields the slope of the function at any given point. So I decided to create a little demonstration, and this is probably one of the first things that got me seriously interested in math.

The black curve is the function (you can edit it from the link above), the red line is the tangent, and the blue line is the normal.

Generalized ellipse functions

Try it yourself:

This graph shows the functions that yield the $x$ and $y$ coordinates of a point as it moves around an ellipse. This is analogous to the $\sin$ and $\cos$ functions for a circle. In fact, if you set $a=c_0$ and $b=c_0$ in the graph above for some $c_0\neq 0$, the resulting functions will be the sine and cosine.

The blue curve is an analog of cosine, and the red is an analog of sine.

I don’t know of any applications for this, but it is pretty interesting to see how the shape of these functions change depending on the characteristics of the ellipse. For example, from the image above, you can tell the major axis is in the $x$ direction because the cosine analog has a greater range than the sine analog. The top and bottom of the blue curve are pinched together because the direction of the point moving along the circle is rapidly changing since the ellipse is wider than it is tall.


Try it yourself:

The blob moves along the function (orange curve) and always stays between the red curves.

Although it has no purpose, this is still one of my favorites. Basically, you can select 1 of 2 blob types. Both blobs follow the given function $f(x)$ as $c$ changes. Press the play button on the $c$ slider or simply move it back and forth manually to see the blob move. The difference between the two blobs is the nature of their shapes. The second blob (“not so blobby blob”) is easier to explain, so I’ll start with that. The general form of a circle centered at $(a, b)$ with radius $r$ is


I wanted the blob to move along the $x$ axis as $c$ varies so I set $a=c$. The blob has to follow the function, though, so I set $b=f(x)$. Note that $f(x)$ is a function of $x$, so chances are it isn’t a constant. When it is, the blob is just a circle. Otherwise, the blob is a sort of ellipsoid shape that points in the direction of the function. Still, this doesn’t feel blobby enough. So began the first blob type (“blobby blob”). I wanted to give the blob an expanding and contracting motion, sort of like a jellyfish. Luckily, math gives us the perfect tool for this: sinusoidal functions. I added in $\cos(x)$ to the $x$ coordinate to create the jellyfish effect and it works pretty well. You can make the blob more interesting by composing and summing additional sinusoids, but it slows down the blob due to greater computation time, which ruins the effect.

Contour Plotter

Try it yourself:

Recently in my math class we went over “contour plots” which show the curves of intersection between a function $f(x,y)$ and a number of constant planes of the form $z=c$. In other words, they show all the points $(x,y)$ such that $f(x,y)=c$. The graphs of functions of two variables are surfaces in three dimensions where the entire $xy$-plane (or a subset of it) is the domain and the output is a single number $z$ which is plotted on the vertical axis. Although Desmos does not support 3D graphs, they do support functions of multiple variables. This allows me to take $f(x,y)-c=0$ where $c$ is a list to make a contour plot. (If you aren’t familiar with how lists work in Desmos, visit

Additional Graphs