Midterm Project: TheThree-Body Problem ENGR 3 Due: 11:59 PM, Sunday, May 3 1 Introduction You have probably learned, perhaps repeatedly, the canonical equation of Newtonian physics: F = ma. This equa- tion underlies essentially every dynamical system, from the interaction between atoms to the movement of entire galaxies. Sometimes, when we’re talking about a particularly simple solution, we can use this simple rule to write out the trajectory of an object. For example, if we throw a ball through the air and ignore air resistance, we know that its height as a function of time can be represented as: y (t) = y0 + v0t+ 1 2 gt2. This is known as a closed-form solution—a solution where we can analytically solve for the position at any time just by looking at the equations. Let’s stick with gravity but up the scale: think about the Earth orbiting the Sun. The acceleration of an object of massm under the influence of gravity from an object of massM is given by: a = GM r2 , (1) where G is a constant and r is the distance between the two objects’ centers. If we consider just the Earth and the Sun, we can again write out a closed-form solution for the Earth’s position at any time in the future! Our solar system is, of course, not just the Sun and the Earth, so let’s also add in the Moon as well. The Moon is quite small compared to these other two bodies, so we expect that it shouldn’t have much impact on the trajectories. And this is true, but it turns out the effect it has is, while small, monumental in what it does to our ability to solve for the solution. That’s because as soon as we add the Moon—our third body—there no longer exists a closed-form solution for the trajectories of any of the three objects. Since we can’t write out an analytic expression for the trajectory anymore, we instead turn to numerical methods, simulating the trajectory one time step at a time to construct the resulting trajectory. We’re going to learn more in the second half of the course about how to utilize some built-in tools in MATLAB to numerically simulate such a system, but for the midterm project, you’re going to create a rudimentary solver and use it to simulate a model of the motion of the Sun, Earth, and Moon. It may sound complicated, but we will tell you step by step what you need to do. 2 Problem Statement You can find on Gauchospace a file called “3body.mat.” Loading the file in MATLAB will provide you with five 1× 3 vectors, corresponding to values for three different objects (so the first value in each vector is for the first object, etc). They are: • masses: the masses of the objects; • x: the initial x-coordinates of the objects; • y: the initial y-coordinates of the objects; • vx: the initial x-velocities of the objects; and • vy: the initial y-velocities of the objects. Using these provided values, you will be responsible for simulating the system. You will do this by writing two functions and a script, the details of which can be found below. The files you will be expected to submit are: 1 • compute_acceleration.m: a function to calculate the gravitational accelerations on each object; • euler_step.m: a function to take a single time step in the simulation; and • threebody.m: a script to load in the initial values and simulate the system, plotting the positions of the three objects at fixed intervals. Finally, Section 4 of this homework asks you to answer some questions. Please submit this as a PDF file generated by Word in addition to your Matlab files. If you don’t have Microsoft Word on your computer, not to worry. It is easy to create a Word document on Box. Instructions on how to do this are included for download in this midterm project on GauchoSpace. 3 Program Specifications Detailed instructions for each component of the midterm project are provided. If you have any questions, don’t hesitate to ask! 3.1 compute_acceleration.m This function should take as input the masses, x-coordinates, and y-coordinates of the three bodies and return two vectors, representing accelerations in the x direction and y direction for each of the 3 masses. The acceleration on an object is equal to the sum of the contributions from the other two objects, according to the gravity equation. Since we need the components of the acceleration in each direction, it is helpful to consider the distances in each direction. Let’s consider the acceleration applied to object 1 with massm1 by object 2, with massm2. If the distance between object 1 (at coordinates (x1, y1)) and object 2 (at coordinates (x2, y2)) in the x-direction is ∆x = x2 − x1, and the distance between them in the y-direction is∆y = y2 − y1, then the total distance between them is: r = √ ∆x2 +∆y2, (2) and the acceleration in each direction is given by: ax = ∆x m2 r3 , (3) ay = ∆y m2 r3 . (4) The total acceleration on object 1 would be this plus the acceleration induced by object 3: ax,1 = (x2 − x1)m2√ (x2 − x1)2 + (y2 − y1)2 3 + (x3 − x1)m3√ (x3 − x1)2 + (y3 − y1)2 3 . (5) Now, using the same reasoning, write down the total acceleration on object 2 and on object 3 and then use all of the accelerations to write your compute_acceleration.m function. To verify your code, we will provide the results for the first acceleration computation: Object 1 2 3 ax 0.000100907029478458 -0.0296000000000000 -0.0672108843537415 ay 0 0 0 3.2 euler_step.m This function takes as input the x and y coordinates, velocities, and accelerations as well as a time step dt, and returns (outputs) the new x and y coordinates and velocities. The new position is calculated by approximating the velocity as a constant over the small time interval dt, and the same is done with the acceleration; this method of numerical 2 simulation is called the forward Euler method. (By the way, Euler is a German name that is pronounced like ‘Oiler’.) The new position is given by: xnew = x+ dt ∗ vx, (6) and the new velocity is given by: vxnew = vx + dt ∗ ax. (7) Using the same reasoning, derive similar equations to update y and vy . Do this for all 3 objects. After the first update using the provided value for dt (see below), the new values should be: Object 1 2 3 x 0 10 10.5000000000000 y 0 0.000273861278752583 0.000344571956871238 vx 5.04535147392290e-06 -0.00148000000000000 -0.00336054421768708 vy 0 0.547722557505166 0.863950323522004 3.3 threebody.m This script will first load in the initial values. This script should also set the value of dt, which you will initially set at dt = 0.0005, and should set t = 0. It should use your function euler_step.m to simulate the 3-body problem from one t to the next, t+ dt. You should continue the simulation until one of two conditions is met: • t reaches 1000, or • the distance between objects 2 and 3 exceeds a value of 3 units. In addition, every 0.5 time units, you should call the provided function “draw_3body” to animate the system. This function takes as inputs the 3× 1 vectors representing the x- and y-coordinates. 4 Analysis Once you have gotten your code working for the provided value of dt, try increasing and decreasing the value of dt. How does changing this affect the simulation? What does this suggest to you about the forward Euler method? Write a brief response to these questions. Later in the course, we will utilize numerical solvers that avoid some of the problems that you have (hopefully!) identified with the forward Euler method. 3
