Homework 5

Due: Tuesday, April 2 at 11:59pm

Submission: On Courseworks

Skeleton Code: hw5.zip

Instructions

Please read the instructions carefully. An incorrectly formatted submission may not be counted.

In this homework, we will start from some skeleton code. You are given the nbody module from the previous homework. We will implement both part 1 and part 2 of this homework on top of the skeleton code, just as in homework 4. To submit, place this module inside a folder called uni-hw5 (for for me this would be tkp2108-hw5). Then compress this folder into a .zip file. For most operating systems, you should be able to right click and have a “compress” option in the context menu. This should create a file called tkp2108-hw5.zip (with your uni). Submit this on courseworks.

Part 1 - Recursion

Recall our discussion on recursion from class. For the following two functions, implement each of them recursively. You may modify the function signatures and/or create helper functions, but be sure to not use global variables.

def exp2(n):
    '''return the value of 2^n implemented recursively. n>=0'''
def isPalindrome(w):
    '''given a string w, return True if it is a palindrome, otherwise return False'''

Here are some example outputs based on my implementations:

from nbody import exp2, isPalindrome

print(exp2(4))  # 16
print(exp2(10)) # 1024
print(isPalindrome("abcba"))  # True
print(isPalindrome("abcb"))  # False

Part 2 - nbody

We’ll be building incrementally on top of the previous homework. You are welcome to start from your version of hw4, or from the solutions. You are also welcome to implement code using pure python, or by using numpy where appropriate. Note that I have added some extra code to generateObjects to build a more realistic universe. It includes initial velocities and a large mass at the center of the system.

We’ll be doing three things in this part:

centerOfMass

First, we want to collect some more info about our data points (to build a more realistic universe model) and in order to annotate our plots with the center of mass of the system. So let’s implement a function to calculate the center of mass.

def centerOfMass(objects):
    """Calculate the center of mass of the objects"""

As a reminder from physics, to calculate the center of mass in a dimension, we sum up the position times mass for each point, then divide by the sum of all the masses.

sumOfForces

Now, we will take our array of objects, and for each other object, sum up the force vectors acting on it from each other object.

def sumOfForces(objects):
    '''Given an Nx5 array of objects of the form (x, y, mass, xvel, yvel), return an Nx2 array of the cumulative force acting on each object in x and y'''

There are many ways to implement this, but don’t overthink it! A for loop inside a for loop will do just fine.

Updating velocity and position based on force

Finally, we will implement the function to update an objects position and velocity based on the forces acting on it over some (frozen) period of time).

def nextPositionAndVelocity(vec, forcevec, delta_t):
    '''Calculate the next velocity and position of vector based on the forces acting on it'''

This code is essentially the same as the trajectory example from the lectures, but broken down into separate x and y steps. So to find the next position and velocity:

# Force = Mass * Acceleration
# Next_V = Current_V + Acceleration * ∆t
# Next_S = Current_S + Current_V * ∆t

# So in x (for example)
next_v_x = current_v_x + f_x / mass * ∆t
next_s_x = current_s_x + current_v_x * ∆t

Skeleton Code

Skeleton code is provided with some function signatures provided, but not implemented. Your job is to implement these functions however you like, but such that the results can be used according to the instructions above. Here is a sample python session from a Jupyter notebook using the library.

Skeleton Code: hw5.zip

import numpy as np
from nbody import centerOfMass, nextPositionAndVelocity, sumOfForces

print(centerOfMass(np.array([[0, 0, 4], [8, 8, 4]])))
print(centerOfMass(np.array([[0, 0, 9], [0, 10, 1]])))
(4.0, 4.0)
(0.0, 1.0)
# Pretend we have an object at position (0, 0)
# with mass 100
# and initial velocity 100 in the x and 0 in the y
obj = [0, 0, 100, 100, 0]


# If we are acting on it for 0.05s
# with a force of -100 in the x
# and 100 in the y
nextPositionAndVelocity(
    obj,
    [-100, 100],
    0.05
)

# the expected output is that in 0.05s
# it moves 5 in the x, 0 in the y
# and has new velocity 99.95 in the x
# and -0.05 in the y
print(obj)
[5.0, 0.0, 100, 99.95, 0.05]
# First, create two fake objects
vec1 = [0, 0, 1935, 0, 0]
vec2 = [3, 4, 1936, 0, 0]

# Then look at their sum of forces,
# they should be equal and opposite
forces = sumOfForces([vec1, vec2])

print(forces[0])
print(forces[1])
[6.000629057280002, 8.000838743040001]
[-6.000629057279998, -8.000838743040001]

Numerical Notes

Because our objects occupy zero space, the skeleton code cheats physics a little bit. This is so we produce prettier pictures in the end. See if you can find it!