Midterm Project: The Three-Body Problem
matlab工程代写 1 Introduction You have probably learned, perhaps repeatedly, the canonical equation of Newtonian physics: F = ma.
You have probably learned, perhaps repeatedly, the canonical equation of Newtonian physics: F = ma. This equation 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:
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 matlab工程代写
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).
- 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:
- 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;
- compute_energy.m : a function to compute the total energy of the system at a particular time
- 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 matlab工程代写
Detailed instructions for each component of the midterm project are provided. If you have any questions, don’t hesitate to ask!
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 mass m1 by object 2, with mass m2. 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:
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 simulation is called the forward Euler method. (By the way, Euler is a German name that is pronounced like ‘Oiler’.)
3.3 threebody.m matlab工程代写
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 t reaches 100.
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.
This function should take as input the masses, x-coordinates, y-coordinates, and velocities of the three bodies and return the total energy of the system.
One way to investigate the accuracy of a numerical method is to check how it treats the energy of the system – we expect the computed solution to satisfy the conservation of energy. For example, we expect that (energy of system at initial time) = (energy of system at final time). If there is a non-neglible difference between these two energy levels, we suspect some numerical error is the cause. Therefore, we are going to develop a function to compute the energy of the system, so that we can use it to check how well our solution method is performing. Let’s consider how to compute the total energy!
The total energy for a body is given by the sum of its kinetic and potential energies (in this case, gravitational potential energy).
Once the energy is calculated for each body, the total energy of the system can be calculated as: Etotal = E1 + E2 + E3 .
Hint: To check if your function is working, you can check the total energy calculation for the initial condition. This should give you Etotal = −0.0045.
4 Analysis matlab工程代写
Once you have gotten your code working for the provided value of dt, try increasing the value of dt and monitor how this change affects the total energy at the final time. For example, try dt = 0.005, dt = 0.05, and dt = 0.5.
Next, take a look at the function “leap_frog_step” , which was provided with the midterm. This function takes the exact same inputs and outputs as your euler_step.m, but performs the Leapfrog method, a slightly different numerical method than the forward Euler method. The only difference is that leap_frog_step also takes a 3 x 1 vector of the masses as an input in addition to the other inputs.
In your threebody.m script, try calling leap_frog_step instead of euler_step.m to run the simulation (comment out the code that calls the euler step, and insert code that calls leapfrog instead). Note: You do not need to modify this function in any way, you only need to call it.Then, investigate the total energy at the final time using this leapfrog method instead – carry out the same process that you did for forward Euler, by trying out different values of dt.
To help with your analysis, fill out the following table with the values of total energy at the final time for both the euler and leapfrog methods, for different values of dt.
Along with filling out the table, write a brief response to the following questions:
1.How does changing dt affect the simulation when using Euler method, and when using Leapfrog method?
2.What does this suggest to you about the forward Euler method?
Later in the course, we will utilize numerical solvers (like leapfrog) that avoid some of the problems that you have (hopefully!) identified with the forward Euler method.
5 Deliverables matlab工程代写
To summarize the deliverables once again, these are the things you need to turn in:
1.compute_acceleration.m: a function to calculate the gravitational accelerations on each object
2.euler_step.m: a function to take a single time step in the simulation;
3.compute_energy.m : a function to compute the total energy of the system at a particular time
4.threebody.m: a script to load in the initial values and simulate the system, plotting the positions of the three objects at fixed intervals.
- This script should call the provided function draw_3body to do the plotting
- Note: when you turn in threebody.m, please turn it in with dt = 0.05, and with the euler_step being called (comment out the leapfrog step again after doing your investigation for section 4).
5. a PDF document with the completed version of the table of final energies for different timesteps and methods, and with a brief written response to the 2 questions in the analysis section