cs222 »

**These are draft notes from subtitles. Please help to improve them. Thank you!**

Contents

- 1 Unit 2: Houston We Have a Solution
- 1.1 01 l Welcome Unit 2
- 1.2 cs222 unit2 01 l Circular Orbit1
- 1.3 cs222 unit2 02 q Geostationary Orbit1
- 1.4 cs222 unit2 02 s Geostationary Orbit1
- 1.5 cs222 unit2 03 p Orbit Error1
- 1.6 cs222 unit2 03 s Orbit Error1
- 1.7 cs222 unit2 04 l Truncation Error1
- 1.8 cs222 unit2 05 q Error Estimate1
- 1.9 cs222 unit2 05 s Error Estimate1
- 1.10 cs222 unit2 06 l Order1
- 1.11 cs222 unit2 07 q Shrink the Error1
- 1.12 cs222 unit2 07 s Shrink the Error1
- 1.13 cs222 unit2 08 l Heuns Method1
- 1.14 cs222 unit2 09 p Programming Heuns1
- 1.15 cs222 unit2 09 s Programming Heuns1
- 1.16 cs222 unit2 10 l Free Euler1
- 1.17 cs222 unit2 11 p Adaptive Step Size CORRECT
- 1.18 cs222 unit2 11 s Adaptive Step Size CORRECT
- 1.19 cs222 unit2 12 l Long Term Behavior1
- 1.20 cs222 unit2 13 l Energy1
- 1.21 cs222 unit2 14 q Speedy Spaceship1
- 1.22 cs222 unit2 14 s Speedy Spaceship1
- 1.23 cs222 unit2 15 p Energy Programming1
- 1.24 cs222 unit2 15 s Energy Programming1
- 1.25 cs222 unit2 16 l Phase Space1
- 1.26 cs222 unit2 17 q Pendular Phase1
- 1.27 cs222 unit2 17 s Pendular Phase1
- 1.28 cs222 unit2 18 l Symplectic Solver1
- 1.29 cs222 unit2 19 p Symplectic Energy
- 1.30 cs222 unit2 19 s Symplectic Energy1
- 1.31 cs222 unit2 20 l Feel The Burn1
- 1.32 cs222 unit2 21 p Boost The Rocket1
- 1.33 cs222 unit2 21 s Boost The Rocket1
- 1.34 cs222 unit2 22 l Save Jim Lovell1
- 1.35 cs222 unit2 23 l Use Of Models1
- 1.36 cs222 unit2 24 q What To Simulate1
- 1.37 cs222 unit2 24 s What To Simulate1
- 1.38 cs222 unit2 25 l Numerics1
- 1.39 cs222 unit2 26 q Numerical Approaches1
- 1.40 cs222 unit2 26 s Numerical Approaches1
- 1.41 28 l Conclusion
- 1.42 cs222 ps2 1 p Symplectic Euler
- 1.43 cs222 ps2 1 s Symplectic Euler CORRECT

Accuracy is an obvious issue with our computer-based solution. To check what the computer does, the easiest way to learn what goes wrong with the computer simulation is to use a toy example, which we can compute by hand as well. This toy example will be a *circular orbit*.

By circular orbit, we mean a spacecraft moving around earth with the earth at the center of the spacecraft's orbit. Let's call the radius of that orbit <math>r</math> and let's call the orbital period the time the ship takes for one orbit, <math>T</math>.

The velocity vector **v** of that spacecraft will be tangent to the orbit. The speed is the length of that velocity vector. We can't easily compute the speed, but we know the speed will be constant. It takes the spacecraft amount of time, capital <math>T</math>, to go around the circumference of the orbit. The distance traveled by the spacecraft is <math>2 πr</math>, the circumference of that orbit. So, the speed <math> s </math> is distance divided by time, or <math> s = 2 πr / T </math>.

The acceleration vector **a** points perpendicular to that orbit, pointing to the center of the earth. The acceleration does not change the speed, that is, the length of the velocity vector. Instead, the acceleration simply changes the direction of the velocity vector.

We can compute the length of that acceleration vector **a** in a similar fashion. This goes as follows. We start with the velocity vector **v** that takes us up and then we move and then the velocity vector changes to the left. At this point, the velocity vector keeps rotating to the left and eventually the velocity vector goes round to full circle.

After one orbital period, it's going to be that blue vector again with phase shift of 90°. So we can't say that the velocity vector executes the same circular motion as does the position vector that we have a phase shift of 90° between them. Nonetheless, we can still apply the same equation. It takes us one orbital period to make one full revolution and so the acceleration vector has a length of <math>2 π*s / T</math>.

Finally, substituting <math>2 πr / T</math> for <math>s</math> from our previous work, the length of the acceleration vector is <math>2π * (2 πr/T) / T</math>, which boils down to <math>(2π)² r/T²</math>. So given the radius <math>r</math> and the orbital period <math>T</math> of this circular motion, we can determine both the speed and the acceleration.

What about the geostationary orbit? The spacecraft should fly around the earth in about 24 hours, so that it seems to be staying above the same spot on the earth. The question is what's the appropriate radius. Give the number of kilometers.

Do you know that the acceleration equals (2π)² r/T², but on the other hand we have Newton's law telling us that this acceleration is G*mE/r². This yields the following. Multiply both sides by r² , (2π)²r/T²=G*mE/r². Let's take r² over here, T² over here, (2π)² over here, so we have, if r³ is whatever remains on the right hand side, r itself is the ³ √ G*mE*T²/(2π)² and if you plug in the numbers, that's about 42,000 km, which is the same as the length of the equator of the earth. So this spacecraft, the distance between the spacecraft and the center of the earth has to be almost the same as the length of the equator of the circumference of the earth. Remember that the distance between that spacecraft and the center of the earth is not the height of that spacecraft. To get the height, you would have to subtract the radius of the earth.

Now that we have an equation for that almost geostationary orbit. Let's look at what the computer simulation does to that, so we start at this point, wait for 24 hours, precisely 24 hours, and look at what's going to happen. Most probably, we're not going to reach that same precise location but end up in a slightly different location, and the objective of this quiz is to compute the distance between the initial position and the position after precisely 24 hours to learn about the error of our method. We've got tons of stuff predefined ready for you to use. Your job is to insert the Euler forward method right at this point. Your job is not to return the positions and the velocities, your job is to return just a single number, namely the error--the distance between the initial position and the final position-- and up front of course compute the step size. Given the number of steps, make sure that you actually are taking, for instance, 200 steps after that initial position or 500 or 1000 and so on.

Of course, the step size simply is the total time divided by the number of steps. To store the positions and the velocities, we have to make space for one more entry, which is the initial condition. We start the position vector at the right location, meaning that its x component is equal to the radius computed above, and the y component of the velocity is equal to the speed computed above. The other components are equal to zero. This is the forward Euler method as we know it. And then we compute the distance between the ultimate position and the initial position using this number function and we return that as the error. So this is the result, what do we learn if we have the step size of 100 seconds, our end result will be off about 1.5*10⁷ m about 15 million meters, and you see that this relationship seems to be linear. We are going to look into that in the next second.

The error that we just saw happening in the quiz is called the global truncation error or GTE if you will. This error stems from using a finite step size. If you would be using an infinitely small step size, you could theoretically achieve zero error. That's part of the truncation using a finite step size and as we saw, this global truncation error depends approximately linearly on the step size. It's more or less some constant times h for the forward Euler's method. This is different for other methods and beware this constant depends on the simulation time. If you look at the global truncation error at this point, it maybe smaller if you take at some further orbital periods, this truncation error may explode. Beware of that this constant is not a constant when it comes to changing the time. When there's a global truncation error, there's of course a local truncation error. That's the error that we're making in every step. See this supplementary material for the local truncation error if you're interested.

So now we know that when we compare the result of the forward Euler method, it's going to deviate from the exact solution. Now assume that you have two numerical solutions. One solution with the forward Euler method of step size h1 and another with a step size of h2 larger than h1, and assume that these two and the exact solution almost line up. The question then is--can we estimate the error between this numerical solution and the exact solution that we do not know based on the difference between those two numerical solutions, which we know. So given the distance between two solutions using the forward Euler method, how would you estimate the error between the exact solution and this first solution with the step size h1? Should it be E times the ratio, E times the ratio of the difference to the first step size, E times the ratio of the first step size to the difference, or E times the inverse ratio?

A sensible estimate would be E times the ratio of the first to the difference. Let's look at why that is the case. We know that for the forward Euler's method, the global truncation error that this kind of error depends on the step size in almost the linear fashion. This error that I am looking for here is the GTE for h1. We see that E, the error between the two numerical solutions, appears here. So now there's two similar triangles. The first one is here and the second one is here. Let's look at the following ratio error/h1. This has to equal E/h2-h1, and hence, the error is equal to E*h1 times up divided by the difference. I'm writing an equal sign here but of course that's just a rough estimate. This dependance is not strictly linear and we're assuming that these three points are almost on a line. This may not be perfect in practice.

Now we can say what the order of a solver is. A solver is a method such as the forward Euler method that we used to solve differential equations on the computer going step by step. For the forward Euler method, we know that the global truncation error, the error after a finite time is something like a constant times the step size. In the general case, it may not be the step size but this step size raise to a certain power, and that power is called the order. The forward Euler method has an order of one, so that we end up with the constant times h. More powerful methods we have larger orders, say two or four, and again, we have to keep in mind that this constant depends on time. It may grow exponentially with time for instance. This constant does not depends on step size, however. Know this site. Be sure not to confuse this notion of order with the notion of order of a differential equation. When we have differential equation of this sort for instance, this is of order two. The order of a differential equation is a different story.

If we had a solver of order 4, remember the forward Euler method is of order 1. If we had a solver of order 4, and if we wanted to shrink the error by a factor of 10, we want to gain one decimal place in accuracy, by which factor do we need to change the number of steps? Should we be using 50% more, 80% more, twice as many, 10 times as many steps?

You should be using 1.8 times the original number of steps. The original error before I change the number of steps would be something like this. Let's call it the old error would be a constant times the old step size to the fourth power. What I want to achieve is that the new error with the new step size is just the 10th part of this, and at the same time, because this is supposed to be a solver of order 4, this is the same constant as before times the new step size to the fourth power. So now we've got two equations that involved the old error and that constant C, this is the first one, this is the second one, we can solve for the new step size which is just basic algebra and you'll find that the new step size is the old step size divided by the ⁴√10 and the ⁴√10 is about 1.8. So we are reducing the step size by a factor of 1.8 which means that we are increasing the number of steps by that very factor of 1.8.

So now we want to actually look at the solver of order 2, the method invented by Karl Heun, Heun's method. It also goes under the name of improved Euler's method. And it is also is one example for a simple Runge-Kutta method in case we haven't got that term. I'll start with the forward Euler method. The position for the first time step is the initial position plus time step times the initial velocity and similarly for the velocity, the velocity is changed by time step times force divided by mass Force divided by mass is acceleration. I have explicitly indicated the dependence of the force on time and position. And what we do now, we do not actually use these results for the next step, we use them as stepping stones if you will. Let's start this result as Xe for Euler and this result is Ve for Euler and the trick now is to use these into immediate result to get a better estimate for the next x and the next v. So how can I improve the estimate for the next position based on Euler's result? And the next position is the initial position plus time step times the average of the initial velocity and the next velocity as predicted by Euler. This takes care of the fact that the velocity is not going to stay constant over the cause of the first time step. Euler's method is based on the assumption that this velocity is almost constant. How can I use this to improve Euler's result? The bad thing about Euler's method is that the velocity actually is not constant. It's not going to stay at the value it has at time equals zero. A better estimate is use the initial velocity and the velocity predicted by Euler and form the average of both. This is what Heun's method does and similarly for the acceleration used to compute the next velocity, I found the average of the initial acceleration F at times zero divided by m and the force after the first time step at the position predicted by the Euler method. For the velocities, I'm using a similar approach. I formed the average of two forces and divided by m to get acceleration. The first force is the initial force times zero, initial position. The second force is the force that I get at the end of the first time step using the position predicted by Euler's method.

So the task now is to actually implement Heun's method, we have provided for you Euler's methods. You can copy these lines and turn them into Heun's method in the end. Again, the error is determined so that we can see by how much Heun's method is better than Euler's method.

I'm working with acceleration here other than the force because acceleration is already offered by the rest of the program, and I said before, I can reuse one of these accelerations, which is why I start the initial acceleration in this variable code, well, initial acceleration. Here comes Euler as we know it, these two lines, and then I'll form the average of the velocities to go to the next step and I'll form the average of the accelerations to go to the next step and that's all. If you run that program, you see two series of data. The upper one is the one that we know, the errors of the forward Euler's method, and the lower one are the errors of Heun's method and you see these are really low. So this method really rocks out. On top of that, you see that this seems to be a parabola. While Euler's method almost forms a line here, Heun's method forms a parabola. So it's good and it even gets better as we shrink the step size. There is a price to pay--we call that acceleration function twice. If you look at what we get for that price, that really pays off.

There is one additional nice thing about Heun's method, it gives us Heun's estimate for the next step and in parallel, we get Euler's estimate for free and we can compare these estimates to find out the error that we are making. By the way, the technical term for this type of method is an embedded Runge-Kutta method if you want to look that up. Remember the two types of errors are introduced before. The one type of error that we're actually interested in, in the end, is the GTE, the global truncation error by how far are we off after a finite's time? What Heun's method allows us is to estimate the LTE, the local truncation error-- the error that we are making at each step. Consider the distance between the Euler result for that step and the Heun result for that step. As we saw earlier, Heun's result is much more accurate than Euler's result. So what we see here is about the error of Euler's method. The LTE of Euler's methods. I'm cheating a little here. I must take the velocity into account. Of course, we have several dimensions, two dimensions concerning x and two dimensions concerning velocity. We have to take all of them into account. So this term as well by how far is the Euler estimate follow velocity of from Heun's estimate? You could simply form the sum here but that would look strange for physicists. These are meters and these are meters per second--that does not work. The most simple idea is to introduce a time in here and time of that simulation. It did makes sense when the velocity is off by a certain amount, but early in that simulation, that error did grow over the course of that simulation. There's other ways of combining these differences--you could form the square here, the square there, and then take the square root but to make things simple, I'm simply using this distance plus the simulation time times that distance. We know that the GTE, the global truncation error for Euler's method is a constant times the step size. We can also show that the local truncation error is another constant times the step size squared. You can check the plausibility of h² appearing here. If I double the step size, I only need half as many steps. Each step has four times the error. I need half as many steps which means that the total error is multiplied by two which we know is the case, but I can do know is to change the step size as I go. I do not need to specify the step size in advance. Once method is able to estimate the perfect step size for a given accuracy, assume that I did the first step with some step size whichever one I chose, but what I want is that this error, the local truncation error of Euler's method is of a fixed size. Let's call that tolerance. What I want is to find the step size. Let's call it h² new, so that this error is off the fixed size called tolerance which would be measured in meters, by the way with this type of error. I can compute C from that very first step. It would be the local truncation error divided by h² old and I can solve for the new step size. So the method actually determines the right step size I should be taking given the tolerance. I can solve for the new step size and find that the new step size is the old step size times the square root of the tolerance divided by this local truncation error that I can compute numerically. So after we did the first step with the old step size, we got an estimate for a perfect step size that produce the local truncation error of such tolerance for Euler's method and we can keep doing that as we go. The step size will be adjusted automatically. We do not have to worry about it. The method would choose the step size for us. One benefit of this method choosing the perfect step size for us, is that we do not have to choose it upfront. The other benefit is that the solution we'll be stepping at appropriate step sizes for instance, when we come close to a planet and the trajectory gets more bent, we need more steps. This method will automatically choose smaller steps as we go and eventually use larger steps again.

Hi! Miriam here making a guest appearance in the middle of the unit. Your task for this quiz is going to be to implement adaptive step sizes. We've already provided Huen's method for you. So you just have to determine the next value of the step size h. There's one important change to the framework, however. You see we're dealing with single vectors here. We're no longer filling arrays of vectors. So the reason that we just are using single vectors is because we don't know upfront how many steps the simulation is going to have. The number of steps is going to depend on the values h that are being computed. So we're not storing a history of x and b anymore because we don't know how many entries that history would have to contain. So when we plot the result, we're not going to be plotting a whole curve. We're going to be plotting just a single dot at each different location. Now you can see that in addition to the step size h right here, we've also created a variable for you called h new. Inside this while loop, we want h to store the current step size and then we'll compute the next step size as h new. At the end of each step, we want to add h to the time it has past, just what's happening right here and I want to set h equal to h new, which signals that we move onto the next step.

What we find here is almost a little interpretation of the equations. However, it's important to notice we compute h new as a modification of h but only use h in the implementation of Heun's method up here. We might run into some difficulties if the numerical error becomes really, really tiny because then as you can see right here in the computation of h new, we're dividing by the error and that would basically make us end up dividing by zero. So one option to prevent this from happening would be to divide not just by the error but by the error plus some really small number. Another precautionary measure that we might want to take would be to limit the range of possible values of the step size. So if we do that, then here's what happens. If h is smaller than 0.1, then the maximum becomes 0.1 and eventually you'll end up with a value of 0.1. Now if time gets larger than 1800, then the minimum becomes 1800 and the new step size becomes 1800 since it limits the range of h new to be between 0.1 and 1800. Now let's look at the output. So we can see that the program takes huge steps along its trajectory and keep in mind that these steps are not at regular time intervals because they represent the steps that we're taking with adaptive step size. So time interval down here is much smaller than the time interval up here. One last thing that we can notice here is that the very first step size is small because we initially set h to 100 seconds and then the method automatically chooses the larger one as you can see right here and then it keeps that larger one, although of course it does change.

We've seen errors appearing on different time scales. The LTE, the local truncation error, that happens at each time step. There is one more time scale that when they look into mainly what happens in the long run, what happens if this system runs for 1,000 years, 10,000 years for the age of the universe, is this solution going to explode or to implode in the long run or will it stick close to the exact solution. The GTE for one final time does not tell us much about that form. One way of looking at the long-term error is for instance to look at energy. If the systems is such that energy should be conserved, we can look at the long-term evolution of energy to learn whether there is an issue.

In this problem, we are concerned with two types of energy, kinetic energy and potential energy. Kinetic energy is due to the speed that is the length of the velocity vector and it grows with speed in a quadratic fashion. It contains a factor of the velocity vector of the spacecraft, length of that vector squared, meaning if increased the velocity by a vector of two, we are going to increase the kinetic energy by a vector of four, which tells us that speeding is really dangerous. There was another factor in here, the mass of the spacecraft and then there is an additional vector of one-half. Potential energy concerns the position of the spacecraft. As we move the spacecraft away from earth, we have to invest energy and that's going to become potential energy. This is the distance of the spacecraft from the center of the earth. When you fire the rockets to get it away from earth, you invest energy and that energy leads to a change in potential energy. In the case of gravity, it's easiest to have potential energy be negative all of the time and only when the spacecraft reaches infinity, it's potential energy go to zero. And the formula for that is pretty much related to Newton's universal law of gravity. We've got the gravitational constant. We've got the product of the masses. Mass of the earth, mass of the spacecraft, and in this case, divided by the first power of the distance. If you remember in the force, we have the second power of the distance. And of course, we have the minus sign in front if we double the distance, we go to half that original value but still with the minus sign and these are the only kinds of energy in our problem and if energy conservation holds, we know the following. The sum of these two energies has to be a constant. It must not depend on time. It has to stay at its initial value and that's something we can use to check our numerical solutions when this expression does not stay yet its original value, we do know that we're in trouble. Note that energy conservation does not hold in all models. In particular, if we have friction and all the heat produced by that friction, energy won't be conserved.

Assume our spacecraft follows this elliptic orbit around the earth. At which point, A, B or C does it have the highest speed?

The highest speed is reached at position B. We have energy conservation in this system that points B the distance of is the smallest, which means that at point B the potential energy is the smallest, But if the sum of potential energy and kinetic energy remains constant, kinetic energy has to be maximum at that point, which means that the spacecraft has achieved its highest speed at this point.

The task now is to study what happens to energy with the forward Euler's method. The method is provided and your task is to fill an array called energy, that's also provided, with the correct values and eventually the results are plotted to see what happens to energy.

To fill the array called energy, we go through all steps which is one more than you may think. The equation as such isn't really surprising. This time for the first time we used a line break, which is a back slash. What the results tells us is that energy keeps increasing. With every orbit, there's a gain in energy and that's pretty much visible in the trajectory, which of course is not the actual physical trajectory. The trajectory keeps expanding and expanding so the long-term behavior of this orbit is pretty poor.

In physics, the conservation of energy is highly related to the conservation of area in phase space. I'm going to introduce that notion in a second. Solvers did also conserve that area, are pretty good at conserving the energy, not exactly but almost exactly. For simplicity, I'm looking at a one dimensional system. The position is one dimension just a single number, the velocity has one dimension, just a single number. Phase space is the coordinate x and v and this is actually use the momentum which is mass time velocity. I don't want you to be confused, so I stick to velocity. For every point in this phase space, you can look it's temporal evolution given the value of x and given the value of v, we can solve the system of differential equations and something is going to happen about that point and another point and to another point. You start from that initial condition to given x and the given v and look at what happens to this x and that v and look at what happens to position and velocity. The trick now is to look at infinitely many initial conditions at the same time. A complete some area in phase space. This is what we're talking about. Let's look at all of these initial conditions. What's going to happen with them? Where is the evolution going? Where is the differential equation taking us? Can this way, we see an evolution of area in phase space? This set of initial conditions will have to say one second end up here and after another second may be end up here and so and so on. An area in phase space remains constant, that's called Liouville's theorem. Why does the area in phase space stay constant as it evolves? Let's look at the first equation. The position changes with the rate of the velocity, so if you rate for a small time and look at wherever ending, x will increase a little here, increase more here, more here, and so on depending on v and this motion does not depend on x, it just depends on v. So, our initial area will be deformed to this which, is the same area. Simply use this part, glue it to that part and you'll see that the blue parallelogram has the same area as the top of the rectangle. So the first equation alone would not change the area. If you look at the second equation on its own, the rate of change of the velocity depends on x and on x only, so if this velocity is changed by that much, also this velocity will change by that much and that velocity will change by that much. The rate of change does not depend on v, it just depends on x, so we could have a different rate of change here but it has to be the same for all velocities and maybe there's a little bit of change here and again it has to be the same for all velocities. So our purple rectangle is transformed into something that may look like this and that still has the same area as the original purple rectangle had. To see that, simply take this slice and push it down and take this slice and push it down and so on and so on. So each equation alone would not change the area. If it continuously apply these rates of change, together we also have no problem. If x evolve for an infinitesimal short amount of time, and then we let v evolve for an infinitesimal short amount of time and then it's the turn of x again and so on and so on--that's going to work. If we apply the first equation for an infinitesimal amount of time, and then we apply the second equation for an infinitesimal amount of time, the area stays the same, and we apply the first equation again and so on and so on. If you step through this and infinitesimally short time steps, that's going to work. Problem is can be applied both equations at the same time and do not use infinitesimally short time steps then it's going to break and that's precisely what the forward Euler's method does. You can see one solution to that problem, you apply the first one for a finite time step then the second one for a finite time step then the first one and second one and so on one after the other, that's going to work, that's going to conserve a phase space area. This leads to the conservation of that area and solvers which have the same property are called symplectic. Sometimes people also use the term geometric integrators. We are going to build the simple most one of these and have a look at the long term behavior and there's something to note all of this works because the force neither depends on time nor on velocity. If the force was influenced by time or velocity that would most probably break this nice behavior. Homework 1 is going to show the evolution of phase space in action producing some nice diagrams.

Let's look at the position and the velocity of the weight that forms the end of a standard pendulum that swings and eventually stop swinging. What happens to a given area in the phase space of that system? Does it remain constant, does it shrinks to half the original size, does it shrink to zero, does it grow to infinity, or does it grow to twice the original size.

The solution is whatever area we're starting this, it's going to shrink to zero. We know that because the motion dies out. Eventually, the position will be fixed and the velocity will be zero. So any initial condition will eventually be taken to the same single point, which means that if I start with any area in phase space and forget all of these initial conditions of a time, they will converged to that single point which means the area becomes zero.

Another one to take the forward Euler's method that you already know so well and turned it into one of these wondrous methods that consist of phase space area. Here I wrote down the equations for standard forward Euler's method. Now, there is a single change that's going to do the job. I do not compute the force at the original position. I compute the force at the new position. This leads to the following. This looks like a programming error on first sight. Soon, I'll be using the initial values on the right-hand side to compute the new values, but actually, it turns out that this the right thing to be doing in contrast to the forward Euler's method. Now, these two steps are executed in sequence. First, we change the position. And you already saw that this does not change the area in phase space and then we take the new position. The velocity remained unchanged and compute the new velocity, which also does not change the area in phase space, which we saw before and then we repeat the first step and so on and so on. This is now a symplectic method. It does not change the area in phase space.

Your job is now to change the standard Euler's method, which we provided here for you, into the symplectic Euler's method and then take a look at the long-term behavior of energy.

There's just a single change in here-- instead of using x at the old step, we're using x at the new step. We see that they're taking almost eight orbits around that planet and the trajectory does not seem to explode or to implode. Even better, we see that the energy follows a periodic pattern. There is no drift up or drift down in energy, and the overall changes of the energy are pretty small, but keep in mind that the energy is not precisely conserve. With some mathematical sophistication, one can show that there's a series of quantities that are close to energy and are conserved better and better.

Now that we know how to reliably solve differential equations, we can return to rockets. What make things easy is that we're going to fire the rocket only for seconds, not for minutes. Think of what would happen if you could give the spacecraft an ultra-short but ultra-strong boost that would lead to an instantaneous change of the velocity. After that of course, also the position changes. We have different trajectory but the only thing that we need to take account of to incorporate that boost is to change the velocity.

Your job is to boost the rocket 2 hours after start by 300 m/s. It maybe a good idea to use this variable that starts whether or not the boost has occurred yet. We provide the rest, namely, Heun's method in this case to provide accurate results.

One way of solving this is the following. We check h times the step number, is the current time. If the current time is larger than or equal to 2 hours, and if we didn't fire the boost before, we fire the boost and then remember that we actually did fire the boost. So next time around, the boost is not going to be fired. You could in principle also check for equality over the time exactly equals to ours or not. But that's very dangerous with floating point numbers. You may miss by a very tiny fraction when the boost is to be done, we change the velocity by 300 m/s that is times the vector which looked somewhat crazy. What we need is a vector of length 1 in the right direction, that is in the direction of travel. 300 times that vector is going to provide us with a vector of length 300 in the direction of travel. So what happens to velocity? We take the existing velocity divided by it's length. This way, we get the vector of length 1 in the direction of travel. 300 times this vector is just the vector that we need to change the velocity vector. There is one more line in this solution. If we fire the rocket, we place a red dot and this is the result. We see the orbit starting here and then the rocket boost takes us to a different orbit.

Now, everything is in place so that you can compute how to rescue Apollo 13, which is homework too. You can use sine and cosine to make the moon go around the Earth. You can use the differential equation to describe the motion of the spacecraft. You know how to fire the rocket to decelerate the spacecraft to take it to a low trajectory around the moon and you can compute by how much to boost the rocket to guarantee a safe return to Earth. The method to be using for that would be Heun's method because we need a quite amount of accuracy. This diagram is really not to scale the Earth and the moon are far more distant as shown here. This figure 8 loop is much more narrow. The packing orbit is pretty close to the surface of the Earth. And this loop around the moon is much closer around the moon as you see here. Actually, for the stimulation, we have to make this loop around the moon a little wider to not spend that much computation time on that.

I want to close with some general comments about the uses of models and simulation. All simulation of Apollo 13 is just one example of what can be done with these techniques from mathematics and computer science. These techniques allow us to do experiments that can't be done in real life. They may be too expensive in real life. For instance, building 1000 cars to find the right aerodynamic shape. They may be too dangerous. For instance, if we ask the astronauts to actually try out different boost or they may be outright impossible--for instance, we want to look at the universe that has different fundamental constants of nature. We may use simulations for simple prediction, for optimization, to study causes and effects for instance of earthquakes or climate change, and we may use models and simulation to determine unobservable quantities such as the composition of the interior of the earth.

What do you think--which of these four questions could be answered with the help of simulation? What's the time of sunset tomorrow, are dolphins more intelligent than elephants, can neutrinos travel faster than light, or who invented the telephone?

What's the time of sunset tomorrow? For sure, we can answer that by simulation. We already simulated the motion of the moon around the earth and in the similar way, we can simulate the motion of the earth around the sun and the rotation of the earth. Are dolphins more intelligent than elephants? To answer this question by simulation which we try our model of how intelligence develops and depends on biology, that is pretty hard to do and may be impossible. Can neutrinos travel faster than light? That's something that physicists have to measure. You could imagine the model that contains more fundamental objects than neutrinos super strings for instance and the neutrinos could be simulated using super strings, but that's again far from the current state of the art and may be impossible. Who invented the telephone? To answer that question, you would have to model all of human history inside the computer That doesn't sound plausible.

And this course we are mostly using the computer to do simulations numerically, that is by number crunching. The major reason for that is that most problems can't be solved by clean equations, unlike school book problems. The real road is way too complex to do that. A second reason is that particularly we have to base our simulations on huge numbers of measurements be it climate's data, be it financial data, be it data from biology, but the numerically approach also has a series of drawbacks. In particular, it's very hard to prove general results using numerics. Is there no way that a particle in some system escapes to infinity? You can look at a single trajectory at a time, but you cannot look at an infinite number of trajectories and numerical methods almost always introduce discretization errors also known as truncation errors. In particular, we have to work with finite steps in time.

For which of these problems that a numerical approach be appropriate: to determine the heat capacity of the oceans, to prove that π/4=1-1/3+1/5-1/7 and so on and so on up to infinity, or to determine whether a control system is unconditionally stable or not.

The numerical approach is just about perfect to determine the heat capacity of the oceans. In particular, as we need tons of geographic data, the numerical approaches are not good to prove something about mathematics, and particularly not when it comes to infinite numbers of operations, and to determine whether a control system is unconditionally stable or not is also not a great task for numerical approaches. We should first do some mathematics to get some general statements about that control system and then we can possibly use numerics on a very abstract level, for instance, to analyze matrices.

Now we come to the first problem of Unit 2. This time we're still dealing with orbits. We have a pendulum and we want to create expressions for the position, velocity, and acceleration. As always, we've given you some code to help you out. Here you can see the time set, the magnitude of the acceleration due to gravity, and the length of the pendulum, which is just the length of the string to which the bulb of the pendulum is attached. What we're asking you to do is to first fill in this definition of the acceleration of the pendulum showing how it depends on the position of the weight. Now if you think about the way that a pendulum swings, you can imagine that if we extend the trajectory, we would get a circle. So position is the length of this curve path. The next thing that you need to do is to fill in this function called symplectic Euler in which logically you will use symplectic Euler method to calculate both distance and velocity. To help you out a little bit, here is a refresher on what the symplectic Euler method says. Another important piece of information if you're not super-comfortable with physics stuff is Newton's second law right here showing the relationship between force, mass, acceleration. So looking back at our code, you can see that we've created empty arrays for you--for position and for velocity. It's up to you to fill these arrays in including the initial conditions. Remember you'll need to initialize x and v where x is zero and v is zero equals something that you're going to figure out. I'm going to give you a few hints though--you can see that we've created this constant called num_initial condition and set it to equal to 50. What you're really doing overall in this problem is looking at 50 different pendulums which each have different initial values for x and v, and to give you a visual of what this looks like, I'm going to show you the final plot that you'll get with this program. So here's the set of graph that you should get as your final result if the program is working correctly. You're going to ignore these top two graphs for now. Let's focus on this bottom plot. This is showing velocity on the vertical axis and position on the horizontal axis. If we look at the top two graphs, we can see that this green set of points right here corresponds to what's happening at time zero. Since green ellipse down here is showing the set of initial conditions for the 50 different pendulums that we're looking at. This green ellipse is kind of like a snap shot of what's happening to all these different pendulums at time zero. The way that I want you to figure out how to set the values for each pendulum for x and for v is just think about this ellipse. Now in Unit 1 you learned about orbits and you'll actually be able to use some of that knowledge in this problem knowing that this green shape is an ellipse. You can see that its major axis right here lies along the line v=0 and its minor axis lies along this line x=2. Now the half length of the major axis is 0.25 so I guess actually it don't looks like the major axis here. It's really the minor axis if these two sets of axes had the same scales applied to them. Either way, this rightmost point on ellipse corresponds to x=2.25 and the leftmost point corresponds to x=1.75. In terms of v, we have values ranging from -2 down here to 2 up here. So think about what equations for x and v you'll need to create an ellipse with these dimensions. That would be how you set the initial conditions for x and v. I think this is a very interesting problem so I hope you enjoy doing it. Good luck.

Let's go over the solution to this problem starting with the definition of the acceleration function. We know the acceleration due to gravity points downward. So let's put this factor into two components. We have one component right here that is parallel to the string of the pendulum and another component that is perpendicular to this green one. Now we know that the acceleration in this direction is going to be exactly cancelled out by the acceleration due to the tension in the rope. So that means that the acceleration we're looking for is really just this pink component, which any point along the path is going to be tangent to the trajectory. Now if we call this angle θ right here, we can figure out the length of this pink component by just saying that it is equal to length of the resultant vector times the sine of θ. The θ down here is actually exactly equal to θ in the diagram of the pendulum itself. So that means to figure out the length of this component, we can use information that we already know about this larger diagram. Since position is just the arc link right here of this imaginary circle, then the measure of that angle in radiant is going to be equal to the length of the arc that it corresponds to divided by the radius of that circle. So that means that in our case, θ is equal to arc length over radius or position over length. So since θ equals position over length and we want the sine of θ, we fill in our definition for acceleration as -g or negative magnitude of the acceleration due to gravity times sine of the position over length. Okay, moving on towards symplectic Euler function. We have to fill in this for loop with the input num_initial conditions. As I said in the InterVideo of the problem, we wanted the initial x to vary from 1.75 to 2.25 and the initial v to vary from -2 to 2 corresponding to the coordinates of every point along that green circle that I'd showed you. Now a convenient way to make a variable cycle through values that are symmetric about an equilibrium value is to use sine or cosine. So we're going to keep that in mind. Now as you learn from the circular orbit problem of Unit 1, if we consider any point along the circumference of the circle, then we can define an angle that corresponds to that point. These are coming from right here as a zero radian mark. You can write the coordinates of this point then as the radius of the circle times the cosine of the angle. That's for the horizontal component and for the vertical component, we get the radius times the sine of the angle. In the phase based plot that I showed you in the InterVideo, we saw that position lying along the horizontal axis and velocity lying along the vertical axis. So we wanted to plot the coordinates of the points on that green circle--the initial condition circle where the position is going to correspond to cosine and the velocity is going to use sine. Now I created a variable called phi. You could pick any name you want I guess. And phi effectively split the circle into 49 segments by marking out 50 different points along the circumference. So every time I increases by one, we're going to step to the next point along the circumference. Since as we saw in the phase base plot, we have a complete circle of green points. The x values of those green points vary like this with 2 as the middle value and the v coordinates vary like that. You noticed that the amplitude in either case corresponds to the half link of the green shape in that direction. So actually we have in a phase base plot is an ellipse for that set of initial conditions. Now that we have our starting additions figured out, we can finally use the symplectic Euler method to proximate the values with x and v at later sets. This code right here is just a direct transition pretty much of the equations that I showed you earlier. Now let's go back to looking at the plot that we get things plugged in but first let's look at our top two plots. The horizontal axis in both of them represents time measured in seconds. The vertical axis in the top one is x measured in meters and here it is v measured in meters per second. So you can see that our initial values of x go from 1.5 to 2.25 and v from -2 to 2. So that corresponds to this green ellipse right here. The most important thing to notice about this bottom graph, which like I said earlier represents phase base is that if we look closely at each one of these ellipses they all have the same area. Now let's look at the shapes that we have down here in this bottom graph. If you look closely and do a bit of calculating, you'll notice that all these different color shapes have the same exact area even though they are well shaped very differently. This is a great example of how phase base is conserve in the system where energy is conserve. Now the fact that its conservation principle holds in this diagram shows how the symplectic Euler method improve upon the accuracy of the forward Euler method. When the forward Euler method is used, it often result in the energy suddenly increasing. So it means that the area of each of these shapes down here will get progressively bigger. Symplectic Euler method, however, confirms much better the equations of motion in physics, It never reflects exactly radical predictions more accurately. Great job with the first problem in Unit 2.