visualization

A couple of weeks ago I made a post about the classical three-body problem, which involves determining the motion of three masses interacting via gravity over time. In that simulation, the three masses were restricted to a plane. Even though we live in three spatial dimensions, a two-dimensional model for celestial orbits isn’t such a bad approximation – the orbits of masses in celestial system are often within a few degrees of the same plane. In the simulation above I’ve built in a third spatial dimension and modeled a system in which the masses do not stay within the same plane. In this simulation, the green and red bodies are 9 and 13 times more massive than the blue body, respectively. The Mathematica code for this one was a bit long, so I’ve posted it here. ...

An orientable surface is a surface on which it’s possible to make a consistent definition of direction. Most surfaces we encounter – like spheres, planes, and tori (doughnut shapes) – are orientable. When visualized in three dimensions, orientable surfaces have two distinct sides. Non-orientable surfaces, on the other hand, have only one side. From Wikipedia, “The essence of one-sidedness is that [an] ant can crawl from one side of the surface to the ‘other’ without going through the surface or flipping over an edge, but simply by crawling far enough.” At any point on a non-orientable surface it’s impossible to uniquely define, for example, the “clockwise” direction. The GIFs above show two examples of non-orientable surfaces: a Klein bottle and a Möbius strip. Mathematica code [Klein Bottle]: xk[u_, v_] := (-2/15)*Cos[u]*(3*Cos[v] - 30*Sin[u] + 90*Cos[u]^4*Sin[u] - 60*Cos[u]^6*Sin[u] + 5*Cos[u]*Cos[v]*Sin[u]) yk[u_, v_] := (-15^(-1))*Sin[u]*(3*Cos[v] - 3*Cos[u]^2*Cos[v] - 48*Cos[u]^4*Cos[v] + 48*Cos[u]^6*Cos[v] - 60*Sin[u] + 5*Cos[u]*Cos[v]*Sin[u] - 5*Cos[u]^3*Cos[v]*Sin[u] - 80*Cos[u]^5*Cos[v]*Sin[u] + 80*Cos[u]^7*Cos[v]*Sin[u]) zk[u_, v_] := (2/15)*(3 + 5*Cos[u]*Sin[u])*Sin[v] kb[u_, v_] := {xk[u, v], yk[u, v], zk[u, v]} Manipulate[ParametricPlot3D[kb[u, v], {u, 0, umax}, {v, 0, 2*Pi}, PlotRange -> {{-1.8, 2}, {0, 4.5}, {-0.75, 0.75}}, Axes -> False, Boxed -> False, PlotStyle -> {Opacity[0.65]}, Mesh -> {20, 11}, MeshStyle -> Directive[Gray, Opacity[0.65], Thickness[0.003]]], {umax, 0.001, Pi}] Mathematica code [Möbius Strip]: xm[u_, v_] := (1 + (v/2)*Cos[u/2])*Cos[u] ym[u_, v_] := (1 + (v/2)*Cos[u/2])*Sin[u] zm[u_, v_] := (v/2)*Sin[u/2] ms[u_, v_] := {xm[u, v], ym[u, v], zm[u, v]} Manipulate[ParametricPlot3D[ ms[u, v], {u, 0, umax}, {v, -1, 1}, PlotRange -> {{-1.1, 1.5}, {-1.5, 1.5}, {-0.5, 0.5}}, PlotStyle -> {Opacity[0.65]}, Axes -> False, Boxed -> False, Mesh -> {20, 5}, MeshStyle -> Directive[Gray, Opacity[0.65], Thickness[0.003]]], {umax, 0.001, 2*Pi}] ...

Given the starting positions, velocities, and masses of three objects interacting via gravity, the classical three-body problem involves determining the motions of the three particles throughout time. What’s cool about the three-body system is that it’s impossible to solve for the motions of the objects exactly. That is, we can’t write down an equation that describes the system. Instead of finding an exact solution, we solve the system numerically, which amounts to finding an accurate approximation. The three-body problem is an example of a chaotic system, meaning that even a slight change in the starting conditions drastically changes the time-evolution of the system. The GIF above shows a planar (i.e., two-dimensional) three-body system. Mathematica code: G = 1; time = 30; mA = 1; xA0 = 0; yA0 = 0; vxA0 = 0; vyA0 = 0; mB = 1; xB0 = 1; yB0 = 0; vxB0 = 0; vyB0 = 0; mC = 1; xC0 = 0; yC0 = 0.8; vxC0 = 0; vyC0 = 0; soln1 = NDSolve[ {mA*Derivative[2][xA][t] == -((G*mA*mB*(xA[t] - xB[t]))/((xA[t] - xB[t])^2 + (yA[t] - yB[t])^2)^(3/2)) - (G*mA*mC*(xA[t] - xC[t]))/((xA[t] - xC[t])^2 + (yA[t] - yC[t])^2)^(3/2), mA*Derivative[2][yA][t] == -((G*mA*mB*(yA[t] - yB[t]))/((xA[t] - xB[t])^2 + (yA[t] - yB[t])^2)^(3/2)) - (G*mA*mC*(yA[t] - yC[t]))/((xA[t] - xC[t])^2 + (yA[t] - yC[t])^2)^(3/2), mB*Derivative[2][xB][t] == -((G*mB*mC*(xB[t] - xC[t]))/((xB[t] - xC[t])^2 + (yB[t] - yC[t])^2)^(3/2)) - (G*mB*mA*(xB[t] - xA[t]))/((xB[t] - xA[t])^2 + (yB[t] - yA[t])^2)^(3/2), mB*Derivative[2][yB][t] == -((G*mB*mC*(yB[t] - yC[t]))/((xB[t] - xC[t])^2 + (yB[t] - yC[t])^2)^(3/2)) - (G*mB*mA*(yB[t] - yA[t]))/((xB[t] - xA[t])^2 + (yB[t] - yA[t])^2)^(3/2), mC*Derivative[2][xC][t] == -((G*mC*mA*(xC[t] - xA[t]))/((xC[t] - xA[t])^2 + (yC[t] - yA[t])^2)^(3/2)) - (G*mC*mB*(xC[t] - xB[t]))/((xC[t] - xB[t])^2 + (yC[t] - yB[t])^2)^(3/2), mC*Derivative[2][yC][t] == -((G*mC*mA*(yC[t] - yA[t]))/((xC[t] - xA[t])^2 + (yC[t] - yA[t])^2)^(3/2)) - (G*mC*mB*(yC[t] - yB[t]))/((xC[t] - xB[t])^2 + (yC[t] - yB[t])^2)^(3/2), xA[0] == xA0, yA[0] == yA0, Derivative[1][xA][0] == vxA0, Derivative[1][yA][0] == vyA0, xB[0] == xB0, yB[0] == yB0, Derivative[1][xB][0] == vxB0, Derivative[1][yB][0] == vyB0, xC[0] == xC0, yC[0] == yC0, Derivative[1][xC][0] == vxC0, Derivative[1][yC][0] == vyC0 }, {xA, yA, xB, yB, xC, yC},{t, 0, time}, MaxSteps -> 100000] x1[t_] := Evaluate[xA[t] /. soln1[[1,1]]] y1[t_] := Evaluate[yA[t] /. soln1[[1,2]]] x2[t_] := Evaluate[xB[t] /. soln1[[1,3]]] y2[t_] := Evaluate[yB[t] /. soln1[[1,4]]] x3[t_] := Evaluate[xC[t] /. soln1[[1,5]]] y3[t_] := Evaluate[yC[t] /. soln1[[1,6]]] Manipulate[Show[ {ParametricPlot[ {{x1[t], y1[t]}, {x2[t], y2[t]}, {x3[t], y3[t]}}, {t, tmax - 0.5, tmax}, Axes -> False, PlotRange -> {{-0.55, 1.45}, {-0.55, 1.08}}, PlotStyle -> {Red, Green, Blue}, GridLines -> {Table[0.25*x + 0.07, {x, -100, 100}], Table[0.25*y + 0.01, {y, -100, 100}]}, GridLinesStyle -> Directive[LightGray]]}, {Graphics[{Opacity[0.7], EdgeForm[Directive[Black]], Red, Disk[{x1[tmax], y1[tmax]}, 0.03], Green, Disk[{x2[tmax], y2[tmax]}, 0.03], Blue, Disk[{x3[tmax], y3[tmax]}, 0.03]}]}, ImageSize -> 600], {tmax, 6.05, 16.05}] ...

A Fourier series is a way to expand a periodic function in terms of sines and cosines. The Fourier series is named after Joseph Fourier, who introduced the series as he solved for a mathematical way to describe how heat transfers in a metal plate. The GIFs above show the 8-term Fourier series approximations of the square wave and the sawtooth wave. Mathematica code: f[t_] := SawtoothWave[t] T = 1; nmax = 18; a0 = (2/T)*Integrate[f[t], {t, -(T/2), T/2}] anlist = Table[(2/T)*Integrate[f[t]*Cos[(2*Pi*n*t)/T], {t, -(T/2), T/2}], {n, 1, nmax}] bnlist = Table[(2/T)*Integrate[f[t]*Sin[(2*Pi*n*t)/T], {t, -(T/2), T/2}], {n, 1, nmax}] fs[t_, nmax_] := a0/2 + Sum[anlist[[n]]*Cos[(2*Pi*n*t)/T] + bnlist[[n]]*Sin[(2*Pi*n*t)/T], {n, 1, nmax}] Manipulate[Column[{Plot[{f[t], fs[t, nmax0]}, {t, -1, 1}, PlotRange -> All, AxesLabel -> {"t", "f(t)"}, PlotStyle -> {{Thick, Black}, {Thick, Red}}, ImageSize -> 700, AspectRatio -> 1/2.8], Row[{"f(t)=", fs[t, nmax0]}]}], {nmax0, 1, nmax, 1}] ...