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