MAS2008 Scientific Computing: Lab 10
3D Plotting and animations

Instructions

Task 1: The dodecahedron

• 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.

Task 2: Solitons

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

%matplotlib widget   
   

• 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.

Task 3: Stacked balls

• 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).