For this homework, the first question is to write an implementation of the so-called trapezoidal method, which is a second order implicit Runge-Kutta method This method is effectively averaging a forward Euler step and a backward Euler step While generally, both backward and forward Euler provide a lot of error in the solution, trapezoidal does a little bit better As some of you may have noticed, it is called trapezoidal method due to its similarities with the trapezoidal integration method you learned in calculus Heres my Matlab implementation of the trapezoidal method As inputs, it takes an initial condition, a step size, and a final time Functionally, this code is very similar to both backward and forward Euler, its just what we do with each step In particular, this line right here This takes the average of a forward Euler step and a backward Euler step Notice that you have to do a forward Euler step to get x to the next time step in order to see the slope there This is very common in implicit methods Notice that while this line of code is just simply the code for forward Euler, I could have just as easily called the forward Euler implementation that I wrote for a previous homework Below the trapezoidal code, for completeness, I provide the simple harmonic oscillator code as well This is also provided in earlier homework solutions Note that k =2, m = 0.5, and beta = 0 are the parameters that well use throughout this homework Part a asks us to find what x and v will be at time 0.5, starting from [-1, -2], using a timestep of h = 0.1, and the given parameters If we do this experiment using the code I just showed, you get x = -1.3811 and v = 0.6211 at time 0.5 For part b, we needed to generate a 500-point trajectory of the same ODE system with a timestep of 0.01, and see which of the following describes the trajectory best This plot is precisely that trajectory, and as you can see, an ellipse is a good description of this For part c, were supposed to keep in mind that there is no damping in the system that is, that friction is equal to 0 and decide whether the trajectories generated with forward and backward Euler in Homework 5.4 are more accurate or less accurate than the trajectory generated with the trapezoidal method The best way to approach this is to simply plot all three on the same plot If we do that, we get this plot All three trajectories started at this point The blue trajectory is backward Euler, and it quickly spirals in to the origin The red trajectory is trapezoidal, and then does the elliptical pattern it should, and the green trajectory is forward Euler, and spirals out very quickly So from this picture, it does seem that the trapezoidal method is providing a more accurate solution, or a better trajectory for these dynamics You can actually formally prove that the trapezoidal method is order h cubed locally, whereas backward and forward Euler are both order h squared So the red really is a better representation of the trajectory For parts d and e, we need to generate a 5,000-point trajectory of the same ODE system and the same initial conditions using a timestep of 0.01, as we did in part b And we want to see qualitatively if theres a difference between this plot and the plot you generated using 500 points Describe that difference, if there is one; and if there is a difference, we need to try to diagnose that difference In this plot, I have plotted both the 500-point trajectory, in red, and the 5,000-point trajectory in blue As you can see, the trajectory has slightly fattened; that is, the trajectory has gotten a little bit wider and takes up a little bit more of phase space Whats happening here is that the numerical error caused in the trapezoidal method, while much, much smaller than either forward or backward Euler, still exists This numerical error is violating the conservation of energy properties of a simple harmonic oscillator with no friction, and since the system is a conservative system, a symplectic integrator should be used, or some type of integrator that conserves energy If youre brave enough to do the suggested problem, and code up RK4, you can see that this problem, where energy is either being gained or lost in the system, persists even for RK4 and RK5 What this should show you is that, while RK4 is the go-to workhorse in dynamical systems, its not appropriate for every dynamical system Even a simple system like the simple harmonic oscillator with no damping can give RK4 grief Just because you have order h to the fourth global error propagation doesnt mean youll get accurate solutions, especially in the cases where special properties need to be preserved like conservation of energy This is something to be aware of, and a reason to know how to code up an ODE solver, and not just use the out-of-the-box method that comes with Mathematica or Matlab or Numpy, for example So the answer to question d is that the trajectory fattened, taking up more space And the answer to question e is that the numerical error is violating the conservation of energy property of a simple harmonic oscillator with beta = 0 Since the system is a conservative system, a symplectic integrator should have been used