MAS2008 Scientific Computing: Lab 10 3D Plotting and animations

Instructions

The following notebooks and videos are relevant for this lab:

Make a 3D plot of a dodecahedron, as follows.
• Use these imports:

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

%matplotlib widget

Instead of %matplotlib widget, you could also try %matplotlib ipympl (after using pip to install ipympl) or %matplotlib qt (after using pip to install pyqt6). With the qt version, your plot will appear in a new window, which may be hidden behind the VS Code window when it is first created.
• Put $\tau_0=(\sqrt{5}+1)/2$ and $\tau_1=(\sqrt{5}-1)/2$. The vertices of the dodecahedron are the 20 points of the form $(\pm 1,\pm 1,\pm 1)$ or $(\pm \tau_0, \pm \tau_1, 0)$ or $(\pm \tau_1, 0, \pm \tau_0)$ or $(0, \pm \tau_0, \pm \tau_1)$. Make a numpy array vertices of shape $(3,20)$ whose columns are these points.
• Write code to calculate the minimum of the distances $\|v_i-v_j\|$ for $0\leq i < j < 20$, where $v_i$ is the $i$th vertex. Call this number edge_length.
• Now make an array called edges of shape (30, 2), whose rows are the pairs [i,j] with $i < j$ such that $\|v_i-v_j\|$ is equal to edge_length. To avoid trouble from numerical inaccuracies, you should actually check whether the distance is less than edge_length plus some small tolerance. You will probably want to construct edges as a list of lists first, and then use np.array(edges,dtype=int) to convert it to an array of integers.
• Put the above code in a function make_dodecahedron() that returns the pair (vertices, edges). Enter this in the online test system.
• Now do ax = plt.figure().add_subplot(projection='3d') to get an object ax that can be used for 3D plotting. Do ax.axis('off') to turn off the axes. Use ax.scatter() to draw all the vertices as blue blobs. Use ax.plot() to draw all the edges as red lines. (You will probably want to run a loop and call ax.plot() once for each edge.)
• You should find that you can click and drag with the left mouse button to rotate the plot, or with the right mouse button to zoom in and out.

The following partial differential equation is called the Kortewegâ€“De Vries equation: $$\frac{\partial\phi}{\partial t} + \frac{\partial^3\phi}{\partial x^3} + 6\phi\frac{\partial\phi}{\partial x} = 0.$$ It is a mathematical model of water waves propagating in a shallow, narrow channel. As the last term involves $\phi$ multiplied by a derivative of $\phi$, this PDE is nonlinear, which usually means that it is impossible to find exact solutions. However, this particular PDE has some special mathematical features which mean that some exact solutions (called solitons) can in fact be found. In this task we will make animated plots of solitons.

Use these imports:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

%matplotlib widget

(Again, you could use ipympl or qt instead of widget.)

Given real numbers $t$ and $x$, we make the following definitions: \begin{align*} q &= \sqrt{2} & p &= \log(3+2q) \\ r &= x-4t & s &= q(x-8t) \\ T &= 32\cosh(2r-p) + 16\cosh(2s-p) + 16 & B &= 4(1+q)\cosh(r)\cosh(s) + (4q-8)e^{r+s} \end{align*} Using these quantities, we can define functions $\phi_i(t,x)$ for $i=0,\dotsc,4$. One can check that these are all solutions of the KdV equation. \begin{align*} \phi_0(t,x) &= 2\cosh(r)^{-2} & \phi_1(t,x) &= 4\cosh(s)^{-2} \\ \phi_2(t,x) &= 2\cosh(r-p)^{-2} & \phi_3(t,x) &= 4\cosh(s-p)^{-2} \\ \phi_4(t,x) &= T/B^2. \end{align*} Write a Python function phi(t, x) that calculates the vector $(\phi_0(t,x),\dotsc,\phi_4(t,x))$. Your code should closely mirror the above definitions; do not try to collapse everything into a single formula. Make sure that your code is vectorized, so that it works correctly when t and x are arrays of the same shape. (It will be even better if your code handles broadcasting correctly when t and x have different shapes, but that is left as an additional challenge for enthusiasts.) Enter your definition in the online test system.

Now plot $\phi_4(-0.8,x)$ for $x$ from $-30$ to $30$, using at least $1000$ points to ensure that the plot is smooth. Experiment with changing the $-0.8$ to various other values, or plotting $\phi_i$ for some other $i$.

Now make an animation of $\phi_4$.
• Set $t$ to be an array of time values between $-3$ and $3$, and set $x$ to be an array of values between $-30$ and $30$.
• Plot $\phi_4(-3, x)$ as before, but instead of just doing ax.plot(...), do line = ax.plot(...)[0] instead. This will give you an object called line which you can use to update the plot. (You should arrange the details so that you get a figure object fig as well as the axis object ax, and you should call ax.plot() instead of plt.plot().)
• Define a function repaint(i). This should calculate the array of values $\phi_4(t_i,x)$, and then call line.set_data(x, y) to update the plot with the new $y$ values. The repaint() function should always return a tuple of objects that have been updated, so in this case it should return the one-element tuple (line,).
• Now call

anim0 = FuncAnimation(fig, repaint,
frames=len(t),
interval=20,
blit=True)

This should display and run the animation.
• Now make a new animation that shows $\phi_0$, $\phi_1$, $\phi_2$ and $\phi_3$. To make the plots visible separately, you should actually plot $\phi_i+5i$ for $i=0,\dotsc,4$. Matplotlib is not very good at setting the y-axis limits in this context, so you should call ax.set_ylim() to set the limits manually. You should label the four plots, using LaTeX to write the labels. (Remember that when a string contains backslashes, it is best to enter it with the prefix r, as in r'$\sin(\theta)$'.) In the initial setup you will need to create four different line objects, which you could assemble into a list called lines. The repaint() function should then update each of these objects and return tuple(lines). The result of the call to FuncAnimation() should be assigned to a new variable, not the same as the one used for the first animation.

What can you say about the similarities and differences between the four functions?
• For large negative $t$, you can see that $\phi_4$ is very close to $\phi_1+\phi_2$. For large positive $t$ we instead see that $\phi_4$ is very close to $\phi_0+\phi_3$ (which is similar to $\phi_1+\phi_2$ but with the peaks in slightly different positions). Experiment with different plots and animations to illustrate this phenomenon.

Make a 3D plot of 20 spheres arranged in a tetrahedron, as follows.
• Use these imports:

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

%matplotlib widget

• The unit sphere can be parametrised as $(x,y,z)=(\cos(u)\sin(v),\sin(u)\sin(v),\cos(v))$ with $0\leq u < 2\pi$ and $0\leq v \leq \pi$. Use this to make a 3D plot of the unit sphere, coloured red. Modify it to plot the sphere of radius $r=0.25$ centred at $(1,1,1)$. You should follow the example of plotting a torus, which was done in the lecture. Note that when doing 3D plots, it is often necessary to call ax.set_xlim(), ax.set_ylim() and ax.set_zlim() (or do something equivalent), because matplotlib will not always choose the limits correctly.
• Now draw 10 spheres of radius $1/2$ arranged in a triangle in the $xy$-plane. Three of the spheres should be centred at $(0,0,0)$, $(1,0,0)$ and $(1/2,\sqrt{3}/2,0)$; from this you can work out where the others should go (draw a diagram on paper if necessary). You should have one sphere for each pair $(i,j)$ with $i,j\geq 0$ and $i+j\leq 3$.
• Now add another three layers, with six, three and one spheres respectively. The bottom left sphere on the first new layer should be centred at $(1/2,\sqrt{3}/6,\sqrt{2/3})$. From this you can work out where the others should go. You should have one sphere for each triple $(i,j,k)$ with $i,j,k\geq 0$ and $i+j+k\leq 3$.
• For the sphere corresponding to the triple $(i,j,k)$, set the colour equal to the triple $((i+1)/4,(j+1)/4,(k+1)/4)$ (which will be interpreted as red, green and blue values).